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Abstract 

Artificial numerical dissipation is an important issue in large Reynolds number compu- 
tations. In such computations, the artificial dissipation inherent in traditional numerical 
schemes can overwhelm the physical dissipation and yield inaccurate results on meshes of 
practical size. In the present work, the space-time conservation element and solution ele- 
ment method is used to construct new and accurate numerical schemes such that artificial 
numerical dissipation will not overwhelm physical dissipation. Specifically, these schemes 
have the property that numerical dissipation vanishes when the physical viscosity goes to 
zero. These new schemes therefore accurately model the physical dissipation even when it 
is extremely small. The method of space-time conservation element and solution element, 
currently under development, is a nontraditional numerical method for solving conservation 
laws. The method is developed on the basis of local and global flux conservation in a space- 
time domain, in which space and time are treated in a unified manner. Explicit solvers for 
model and fluid dynamic conservation laws have previously been investigated. In this paper, 
we introduce a new concept in the design of implicit schemes, and use it to construct two 
hi ghl y accurate solvers for a convection-diffusion equation. The two schemes become iden- 
tical in the pine convection case, and in the pure diffusion case. The implicit schemes are 
applicable over the whole Reynolds number range, from purely diffusive equations to purely 
inviscid (convective) equations. The stability and consistency of the schemes are analysed, 
and some numerical results are presented. It is shown that, in the inviscid case, the new 
schemes become explicit and their amplification factors are identical to those of the Leapfrog 
scheme. On the other hand, in the pure diffusion case, their principal amplification factor 
becomes the amplification factor of the Crank-Nicolson scheme. We also construct an ex- 
plicit solver with the treatment of diffusion being based on that in the implicit solvers. The 
explicit solver has only a CFL stability limitation on the Courant number, yet it retains the 
second-order spatial accuracy of the implicit schemes. 
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1 Introduction 


The method of space-time conservation element and solution element (the CE/SE method, 
for short) is a new numerical discretization method for solving conservation laws ([l]-[27]). 
The distinguishing features of the method are (i) an emphasis on solving a general space-time 
integral form of the conservation laws, (ii) the requirement of local and global conservation of 
space-time fluxes as a basis for the method, (iii) the construction of numerical schemes so as 
to reflect the initial- value/boundary-value nature and invariance properties of the governing 
equations, (iv) use of only the simplest discretization stencils, (v) the absence of upwind- 
biasing of fluxes and approximate Riemann solvers with their attendant complexities, (vi) 
the absence of directional splitting techniques and the difficulties associated with them, for 
flow computations in multiple spatial dimensions, and (vi) the treatment of flow property 
gradients as additional unknowns, thus eliminating the need for gradient reconstruction 
using ad hoc polynomial fitting and flux limiters. Reference [11] is a review of the CE/SE 
method that details the differences between the CE/SE method and the traditional numerical 
methods for conservation laws. 


Based upon the physics of a scalar convection equation, [6] contains a discussion of the 
requirements for a numerical scheme to be an ideal solver for such a conservation law. The 
desirability of the distinguishing features mentioned in the previous paragraph follows from 
the requirements listed in [6], and from an emphasis in the CE/SE method on simplicity, 
generality and accuracy. With these requirements and this emphasis as guiding principles, 
several two-level explicit schemes were constructed in [6, 4] to solve (i) the pure convection 


equation 


du du 

ai +a ai = 0 


(i.i) 


with given initial values, and (ii) the convection-diffusion equation 


du du d 2 u 
dt +a di~^d^ = ° 


with given initial and boundary values. In the foregoing equations, the convection velocity a, 
and the viscosity coefficient [i (> 0) are constants. These schemes were then extended to solve 
the 1-D time-dependent Euler and Navier-Stokes equations of a perfect gas [6, 4]. Moreover, 
except for the Navier-Stokes solver, the above 1-D schemes have been generalized to their 
2-D counterparts [7, 5]. Because of the inherent simplicity and generality of the current 
method, the above multidimensional generalization is a straightforward matter. Also, as a 
result of the similarity in their designs, each of the above 2-D schemes shares with its 1-D 
version virtually the same fundamental characteristics. 

We now discuss the appropriateness of using an implicit, rather than explicit, solver for 
initial- value/ boundary- value problems such as that represented by Eq. (1.2). For a two-level 
explicit scheme, the value of a solution at any mesh point has a finite domain of dependence 
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at the previous time level. As an example, consider a finite-difference solver for Eq. (1.1). 
Let u", the mesh value of u at any mesh point (j, n) (point P in Fig. 1(a)) be determined 
by Ujli, u " -1 , and Uj+l • Then the numerical domain of dependence of u” at the (n — l)th 
time level contains three mesh points. Also, one can see that ti" is dependent only on the 
initial data given on the line segment AB. 

For an initial-value problem, such as a time-dependent Euler problem or a problem in- 
volving Eq. (1.1), the exact solution at any point in space-time also has a finite domain of 
dependence on the initial plane. As a result, explicit schemes could be ideal solvers for such 
a problem if they satisfy the requirement that the physical domain of dependence be a subset 
of the numerical domain of dependence. The latter requirement manifests itself as a CFL 
stability limit on the Courant number. 

On the other hand, the solution of an initial-value/boundary-value problem at any point 
in space-time is dependent on the initial data and the boundary data up to the time of the 
point under consideration. As an example, consider a problem involving Eq. (1.2). As shown 
in Fig. 1(b), the exact solution at point P is dependent on the initial/boundary data given 
on EC, CD, and DF, where E, P, and F are at the same time level. Let this problem 
be solved using the explicit scheme that was explained with reference to Fig. 1(a). Let P 
also be a mesh point ( j , n). Then the numerical solution value it" is dependent only on the 
initial/boundary data given on AC, CD and DB. It is completely independent of those 
data given on AE and BF. Contrarily, if the same problem is solved using an implicit 
scheme, then u" is dependent on the initial/boundary data given on EC, CD, and DF. In 
other words, the numerical domain of dependence of the implicit scheme is consistent with 
the physical domain of dependence of the problem under consideration, while that of the 
explicit scheme is not. 

Two observations can be made as a result of the above discussions. 

(a) Generally speaking, an explicit scheme is not an ideal solver for an initial- value/ boundary- 
value problem. Because a time-dependent Navier-Stokes problem is such a problem, the 
above argument implies that an explicit scheme cannot be used to solve a time-dependent 
Navier-Stokes problem except for the special circumstance in which errors caused by ne- 
glecting certain initial/boundary data (such as those given on AE and BF in Fig. 1(a)) 
are relatively small. The factors that help achieve the above special circumstance include: 
(i) a small time-step size to spatial-mesh interval ratio, (ii) a small time rate of change 
of boundary data, and (iii) a small contribution of the viscous terms in the Navier-Stokes 
equations relative to that of the inertial terms. Note that condition (iii) may be met by a 
high-Reynolds-number flow. Under the special circumstance when an explicit solver may be 
used, such a solver may have the advantage of being computationally less expensive than an 
implicit solver for the same problem. 

(b) Generally, an implicit scheme is not an ideal solver for an initial-value problem. This 
is because the domain of dependence in the former is far greater than in the latter. As a 
result, although an implicit solver meets stability criteria, the resulting solution tends to be 
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contaminated by extraneous information. 

With these considerations in mind, in this paper we will construct and analyse two implicit 
solvers for Eq. (1.2). For reasons (including remark (b) above) that will be explained in 
Section 3, the implicit solvers will be constructed in such a way that when the viscosity 
coefficient p vanishes, the schemes will become an explicit scheme for Eq. (1.1), termed the 
a scheme, that was described in [6]. A special case of the implicit schemes was described 
without analysis in [8]. We will also construct and analyse an explicit scheme which also, 
in the inviscid case, reduces to the a scheme. The explicit scheme is suitable as a solver 
for Eq. (1.2) either when the time-step to spatial-mesh-interval ratio is small enough, or 
when the time rate of change of boundary data is zero and only a ‘steady-state’ limit of the 
solution is of interest. 

The CE/SE explicit solvers referred to earlier, and the implicit and explicit solvers to be 
described herein, are all built on a foundation of flux conservation in space-time. Although 
the differential equation representing a conservation law is used in the CE/SE method, 
emphasis is placed on a space-time integral form of the conservation law. In the following 
paragraphs, we explain the space-time integral formulation of conservation laws. 

Let Eqs. (1.1) and (1.2) be in dimensionless form. Let xi—x and x 2 = t be considered as 
the coordinates of a two-dimensional Euclidean space E 2 . Consider the integral conservation 
law 

i h-ds = 0. (1.3) 

JS(R) 

Here (i) S(R) is the boundary of an arbitrary space-time region R in E 2 (see Fig. 2), (ii) 
h is a current density vector in E 2 , and (iii) ds = dan with da and n, respectively, being 
the area and the outward unit normal of a surface element on S(R). Note that (i) h ■ ds is 
the space-time flux of h leaving the region R through the surface element ds, and (ii) all 
mathematical operations can be carried out as though E 2 were an ordinary two-dimensional 
Euclidean space. 

If we set 

h = (au,u ) (1.4) 

in Eq. (1.3), then it can be shown that Eq. (1.1) is the differential form of Eq. (1.3), under 
suitable smoothness assumptions on u. This can be done by using Gauss’ divergence theorem 
in the space-time E 2 . If, instead, we set 



du 

{au - /i— ,u) 


(1.5) 


in Eq. (1.3), then it can similarly be shown that Eq. (1.2) is the differential form of Eq. 
(1.3), again under suitable smoothness assumptions on u. In both cases, Eq. (1.3) is the 
more fundamental and more general expression of the physical conservation law. 

At this juncture, note that the conservation law Eq. (1.3) appears in a form in which 
space and time are unified and treated on the same footing. This unity of space and time 
is a key characteristic that distinguishes the current method from most of the traditional 
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methods. It provides flexibility in the shape of the discrete elements over which conservation 
is enforced. Also, it is the conservation of space-time flux that is enforced, rather than 
the traditional conservation of spatial flux with time evolution possibly being outside the 
conservation framework. 

As mentioned above, the implicit and explicit solvers developed in this paper contain the 
explicit a scheme for the pure convection equation as a special case. An understanding of 
the latter scheme is the best route to understanding the properties of the new implicit and 
explicit solvers. Hence, we review the a scheme in the next section. In Section 3, we develop 
the discrete conservation equations for the implicit solvers. The solution procedure for each 
implicit solver is presented in Section 4. The stability of the implicit schemes is discussed in 
Section 5. The consistency and truncation errors of the implicit schemes are investigated in 
Section 6. A modification of the initial conditions to improve solution accuracy is described 
in Section 7. We develop and analyse the new explicit scheme in Section 8. Some numerical 
examples demonstrating the properties of the new schemes are displayed in Section 9, and 
we end with a summary of the paper. 


2 Review of an Explicit Scheme 


In this section, we review the a scheme, which is the inviscid version of the explicit a-p 
scheme [6, 4]. The a scheme is a nondissipative solver for the pure convection equation, 
Eq. (1.1). To achieve consistency of the notations used in this section and the sections that 
follow, the notations used here will be slightly different from those used in [6, 4]. 

Let Qi denote the set of mesh points (j, n ) in E 2 (dots in Fig. 3) where n — 0, ±1, ±2, ±3, 
and, for each n, j = n,n ± 2 , n ± 4 , , . . .. Note that only a finite portion of this unbounded 
space-time mesh is seen in Fig. 3. We select a numerical representation of the solution that 
is piecewise Unear in the ( x , t) coordinates. There is a solution element (SE) associated with 
each ( j, n, ) G Let the solution element SE(j,n) be the space-time region bounded by 
the dashed curve depicted in Fig. 4. It includes a horizontal line segment, a vertical line 
segment, and their immediate neighborhood. The exact size and shape of the neighborhood 
is immaterial to our purpose. 

For any ( x , t) G SE (;, n), u(x, t), and h(x. t), respectively, are approximated by it*(x, t ; j, n ) 
and h*(x,t',j, n) which we shall define shortly. Let 


u"(x,t]j,n) = f u“ + (ti,)”(i - Xj ) + (u, )”(t- t"), 


(2.1) 


where (i) it", (it*)", and (u t )” are constants in SE(j, n), and (ii) (xj,t n ) are the coordinates 
of the mesh point ( j , n). Note that 


u m (x j ,t n ;j,n) = it"; 


du*(x,f ; j,n) 

dx 


= K)”; 


du*{x t t-,j,n) 

dt 


= (u t )]. (2.2) 


Moreover, if we identify it”, (it*)", and (u t )”, respectively, with the values of u, du/dx, 
and du/dt at (xj,t n ), the expression on the right side of Eq. (2.1) becomes the first-order 
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Taylor’s expansion of u(x,t ) at (xj,t n ). Thus, it", (u x )”, and (u t )” are respectively the 
numerical analogues of the values of u, du/dx, and du/dt at ( Xj,t n ). We mention here that 
the choice of basis functions for the linear approximation is both convenient and meaningful, 
because the conservative variable u and its gradient arise naturally in the statement of 
the conservation law. Numerical solution of an initial value problem for Eq. (1.1) then 
reduces to determining the three unknown constants in Eq. (2.1) associated with each mesh 
point (j, n). This will be done by specifying that the numerical approximation satisfy three 
equations associated with each mesh point. These three equations each model the governing 
conservation law. 

To begin with, we shall require that u = u*(x,t]j,n ) satisfy the differential form of the 
conservation law, Eq. (1.1), within SE(j,n), i.e., 

(u t )^ = -a(u x )] (2.3) 

We note in passing that, since u* is a Unear approximation, Eq. (2.3) can alternatively be 
arrived at by requiring u* to satisfy the integral form of the conservation law, Eq. (1.3), 
within SE(j,n). Combining Eqs. (2.1) and (2.3), one has 

u*{x,t;j,n) = u" + (u x )][(x-x j )-a{t-t n )\ , (x,t) e SE{j,n). (2.4) 

From Eq. (1.4), h = (era. u). We therefore define 

h*(x,t;j,n) = ( au*{x,t;j,n ), u*(x,t ; j,n)) . (2.5) 

Let i ?2 be divided into nonoverlapping rectangular regions (see Fig. 3) referred to as 
conservation elements (CEs). As depicted in Figs. 5(a) and 5(b), the CE with its top- 
right (top-left) vertex being the mesh point ( j , n) e ffy is denoted by CE_(j, n) (CE+(j, n)). 
Obviously the boundary of CE_(j, n) (CE + (j, n)) is formed by subsets of SE(j, n) and SE(j - 
1, n- 1) (SE (j -+l,n — l)). The numerical approximation u * is now determined by requiring it 
to satisfy the integral conservation law applied to the CEs. This requirement yields the other 
two equations associated with each mesh point. The current approximation of Eq. (1.3) is 

F±(j,n) = <fi h'-ds = 0, (2.6) 

JS(CE±(j,n)) 

for all (j, n) G flj. In other words, the total flux leaving the boimdary of any CE is zero. 
Note that the flux at any interface separating two neighboring CEs is calculated using the 
information from a single SE. As an example, the interface AD depicted in Figs. 5(a) and 
5(b) is a subset of SE(j, n). Thus the flux at this interface is calculated using the information 
associated with SE(j>, n). 

Because (i) The CEs associated with Qi can fill any space-time region, and (ii) the sur- 
face integration across any interface separating two neighboring CEs is evaluated using the 
information from a single SE, the local conservation condition Eq. (2.6) leads to a global 
conservation relation, i.e., the total flux leaving the boundary of any space-time region that 
is the union of any combination of CEs will also vanish. 
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(2.7) 


With the aid of Eqs. {2.4)-(2.6), it can be shown that (see [6, 4]) 

F±{j, n)/Ax = ±(1 - v 2 ) [«)” + «)£i] + (1 T v) («? - u£l) , 

where v = aAt/Ax is the Courant number, and (u+ )” = (Ax/2)(u I )". Note that here 
A t and At, respectively, represent the same mesh interval and time-step size which were 
denoted by Ax/2 and At/2 in [6, 4]. Eq. (2.6) or Eq. (2.7) represents two equations, one 
corresponding to the upper signs, and the other to the lower. This convention will also be 
used at other places in this report. Using Eqs. (2.6) and (2.7), u" and (u+)", which are 
considered as independent unknowns at the mesh point (j, n), can be solved for in terms of 
u™±\ and (O^ 1 if 1 - z/ 2 7^ 0. It is shown in [6, 4] that, for all (j,n) e fii, 

q{j,n) = Q+q{j -l,n-l) + Q-q{j + l,n-l), (1 - v 2 ^ 0). (2.8) 

Here (i) q{j,n ) is the column matrix formed by u" and (u+)”, and (ii) 


Q± = (1/2) 


1 ± v ±(1 - v 2 ) \ 
Tl -1 ±v )■ 


(2.9) 


Eq. (2.8) defines a marching scheme. Because this scheme models Eq. (1.1) which is charac- 
terized by the parameter a, it is referred to as the a scheme. It is also shown in [6, 4] that 
the a scheme is stable provided that v 2 < 1, i.e., provided the CFL criterion is satisfied. 
Given tt” and (u+)" at each point on a segment of some initial line of constant n, Eq. (2.8) 
may be used to find the numerical solution at points (j, n ) , at a later time. Because of the 
stability limit, such points must lie within the region of influence of the initial line segment. 
T his stability requirement is reflective of the character of Eq. (1.1). 

Because the local condition Eq. (2.3) explicitly allows us to easily eliminate (u t )” in favor 
of (tti)", the other flux conditions (2.7) are written in terms of u” and (u+)”. Thus (u t )™ 
does not need to be stored, and u ” and (tt+)" can be considered to be the basic marching 
variables that are marched via Eq. (2.8). 

The a scheme is in many ways an ideal solver for Eq. (1.1). The property of the a scheme 
of greatest significance for the construction of an implicit solver for Eq. (1.2) is its non- 
dissipative nature. The a scheme is the only two-level explicit solver of Eq. (1.1) known 
to the authors to be neutrally stable, i.e., free from numerical dissipation. It also has the 
s im plest, stencil, i.e., a triangle with a vertex at the given time level and the other two vertices 
at the previous time level. Because the flux at an interface separating two neighboring CEs 
is evaluated using information of a single SE, no interpolation or extrapolation is required. 
Moreover, the a scheme is a two-way marching scheme, i.e., a backward marching scheme in 
which q(j,n) is determined in terms of q(j—l,n+l) and q(j + l,n + l) can also be derived 
from Eq. (2.7). This is entirely in keeping with the invariance under time-reversal displayed 
by Eq. (1.1). The backward marching scheme may be used to find the numerical solution at 
points (j, n), at an earlier time. Such points must lie within the region of influence of a given 
initial line segment. These and other nontraditional features of the a scheme are discussed 
in depth in [6, 4]. 
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In the above construction of the a scheme, we use the SEs and CEs of the mesh points 
marked by dots in Fig. 3 and Fig. 6. A similar construction can be performed by using the 
mesh points marked by triangles in Fig. 6. Let Q 2 denote the set of mesh points (j, n) in E 2 
(triangles in Fig. 6) where n = 0, ±1,±2,±3, . . ., and, for each n, j = n ± l,n±3,n±5,.. .. 
Let the SEs and CEs of fl 2 be defined by using Figs. 4 and 5(a)-(b), with dots replaced 
by triangles. Obviously (i) the CEs of Q 2 also fill any space-time region, and (ii) the a 
scheme can also be constructed using the SEs and CEs of fi 2 . This new scheme is defined 
by Eq. (2.8) with (j,n) € fi 2 . 

We remark on what may appear paradoxical about the overlapping of SEs and CEs 
associated with fh and Q 2 . Let (j,n) € Oi- Then (i) (j + l,n) € fl 2 , (ii) SE(j, n) and 
SE(j + l,n) overlap each other, and (iii) as depicted in Fig. 7, CE+Jj, n) and CE_ (j 4- l,n) 
represent the same rectangle in E 2 . However, because the function h* used in the evaluation 
of F + (j, n) is tied to a pair of SEs associated with Oi, while that used in the evaluation of 
F-(j- 1-1, n) is tied to another pair of SEs associated with fi 2 , F+(j, n) = 0 and F-(j+ 1, n) = 0 
represent two completely independent dux conservation conditions. 

To prepare for the development of the implicit solvers to be described in the next section, 
we shall combine the above two independent schemes into a single scheme. The new scheme, 
referred to as the a(2) scheme, is defined by Eq. (2.8) with (j, n) 6 f2 where Q = f Qi U fi 2 . 
Obviously, a solution of the a(2) scheme is formed by two decoupled solutions with each 
being associated with a mesh that is staggered in time. Several classical schemes also have 
this property. Among them are the Leapfrog, the DuFort-Frankel, and the Lax schemes [28]. 
We mention that the meshes fix and fi 2 are duals of each other, in the sense used in the 
literature of fini te- vol um e methods or of finite-element methods. We term their union the 
dual space-time mesh. 


3 The implicit schemes 

3.1 The a-fi{Il) scheme 

An im plicit, solver for Eq. (1.2), referred to as the a-p(Il) scheme, will be discussed in this 
sub-section. Here “7” stands for “implicit”, and “1” is the identification number. This 
solver is the model for implicit time-dependent Navier-Stokes solvers under development. It 
is constructed to meet two requirements given in the following discussion: 

(a) With a few exceptions, numerical dissipation generally appears in a numerical solution 
of a time-marc hing problem. In other words, the numerical solution dissipates faster than 
the corresponding physical solution. For a nearly inviscid problem, e.g., flow at a large 
Reynolds number, this could be a serious difficulty because numerical dissipation may over- 
whelm physical dissipation and cause a complete distortion of solutions. To avoid such a 
difficulty, the model solver is required to have the property that the numerical dissipation 
shall approach zero as the physical dissipation approaches zero. 
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(b) The convection term and the diffusion term in Eq. (1.2) involve the spatial derivatives 
of first order and second order, respectively. Thus, in a spatial region where a solution is 
very smooth, the diffusion term is negligible when compared with the convection term. As 
a result, the effective physical domain of dependence is more or less dictated by Eq. (1.1)- 
To prevent excessive contamination of the solution by extraneous information, the implicit 
solver shall be required to become an explicit solver in the limiting case in which the diffusion 
term vanishes. 

On the basis of the remarks (a) and (b), an implicit solver for Eq. (1.2) should reduce to 
an explicit nondissipative solver when p = 0. Because of these requirements, the implicit 
solver will be constructed such that it reduces to the a( 2) scheme if p = 0. The former differs 
from the latter only in the extra modeling involving the diffusion-related terms. Note that 
the presence of viscosity is felt through (i) the diffusion term -pd 2 u/dx 2 in Eq. (1.2), and 
(ii) the spatial diffusion flux component —pdu/dx in the flux vector in Eq. (1.5). 

We now construct the a-p(ll) scheme, to numerically solve initial/boundary value prob- 
lems involving Eq. (1.2). Let the initial and boundary values be given as 

u(x, 0) = ti/(x) for 0 < x < JAx, (3.1) 

u(0,t)=u L (t) for t > 0, (3.2) 

u(JAx,t) = u R (t) for t > 0, (3.3) 

where uj(x ), ui(t), and u R (t) are given functions, and J and Ax are defined next. Consider 
the mesh depicted in Fig. 6 (J > 4), which is suitable for problems with stationary bound- 
aries. There are J uniform spatial mesh intervals of size Ax. The temporal mesh intervals 
axe each of size At. The SEs and CEs are defined in the interior as they are for the a(2) 
scheme, which was described in the previous section. It follows that, for the current case, 
(i) and fl 2 are restricted by the conditions n > 0 and J > j > 0, (ii) CE ±(j,n) are not 
defined if n = 0, (iii) CE_ (j, n) is not defined if j — 0, and (iv) CE+(j, n) is not defined 
if j = J . Items (iii) and (iv) imply that only one conservation condition is associated with 
a boundary mesh point. Boundary values will supply the additional information needed to 
determine all discrete unknowns at a boundary mesh point. Obviously, the definition of 
SE(j, n ) also needs to be appropriately modified if j = 0, or j = «/, or n = 0. 

We seek to determine an approximation u* to the solution of the initial/boundary value 
problem. As in the a(2) scheme, we specify that the approximant have a piecewise linear 
variation with space and time, i.e., Eq. (2.1) will still be assumed: 

u*(x,t-J,n)=u] + (u x )^{x-x j ) + {u t )]{t-t n ) , (x,t)€ SE(j,n). 

Numerical solution of the problem then reduces to determining the three unknown constants 
in Eq. (2.1) associated with each mesh point (j,n). At each boundary mesh point (0,n) 
or ( J,n ), the applicable boundary condition will be used to determine two of the three un- 
knowns. At each mesh point ( j, 0) on the initial line, the initial condition will be used to 
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determine two of the three unknowns. Each of the corner points (0, 0) and ( J, 0) is simulta- 
neously a boundary point and a point on the initial line. At each of these corner points, all 
three constants in Eq. (2.1) are determined from the initial and boundary conditions. 

Three equations, each representing a balance of space-time fluxes, will be developed at 
each interior mesh point (j, n). As in the a scheme, these equations will be numerical 
analogues of (a) the differential form of the conservation law evaluated at (j,n), (b) the 
integral form for CE+(j,n), and (c) the integral form for CE_ (j, n). Also, one of (a), (b) 
and (c) will be written at each boundary or initial point which is not a corner point, thus 
providing a third condition at each such point in addition to the two conditions coming from 
initial or boundary values. Thus, the number of equations match the number of unknowns, 
with three equations and three unknowns being associated with each mesh point. 

We proceed to develop the numerical analogue of the differential form of the conservation 
law at the point (j, n). We require the behavior of u* to model Eq. (1.2). Simply substituting 
u = u* in Eq. (1.2), with u* given by Eq. (2.1), would lead to Eq. (2.3) of the inviscid a 
scheme. Because the expansion in (2.1) is only linear, the second-derivative diffusion term 
in Eq. (1.2) is not modeled in Eq. (2.3). Instead, in the present scheme, we replace Eq. (2.3) 


with 


0 u t)j — a { u x)j + 9 A |( W x)j+l ( w x)j-l i J > j > 0 , 


2 Ax L 


(3.4) 


for all n > 0. Eq. (3.4) is the numerical analogue of the differential condition Eq. (1.2). A 
comparison between Eqs. (1.2) and (3.4) reveals that the diffusion term —pd 2 u/dx 2 at an 
interior mesh point G fii (Q 2 ) is modeled by a central-difference approximation involving 
the values of u x at two neighboring mesh points G O 2 (fii). Note that Eqs. (2.1) and (3.4) 
imply that, for J > j > 0 and (x, t) G SE(j. n), 


u*(x,t;j,n) 

= Kj + («,)“(! - x,) + [(u,)“ +1 - («.)*-i] - »(».)■} (t - «"). (3.5) 


Next, we develop the two remaining flux balance equations by requiring the numerical 
space-time fluxes to satisfy the integral form of the conservation law for the CEs defined 
in the previous section. To do this, first we form the numerical analogue of the flux vector 
defined by Eq. (1.5). Eq. (2.5) is replaced by 


h ( X)t',j,Ti ) — (qu (Xjtjj^n) pu x (x,t,j,n), u (x,t,j,n)), 


(3.6) 


for (x, t) G SE (j,n). Here, the function u* x in the diffusive flux is defined as follows. For 
(x, t) e SE (j,n), 


u*0M; j,n) 


def 


w(u x y + w(u x )” 1 , if t < t n ; 
«/(«,)? + u/K)” +1 , if t > t n . 


(3.7) 


where tu and u/ are weighting factors to be specified. The notations tu = \—w and w' = 1—w' 
have been used for compactness. Note that if w' = w, then u*(x, t ,j, n — 1) = u*(x, t ; j , n) 
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in the region of overlap of SE (j, n) and SE (j, n — 1), i.e., in this case u* is uniquely defined 
at any (x, t) where it is defined. 

With the above modifications, the a-p(Il) scheme is defined by assuming Eq. (2.6). Note 
that: 

(a) At a mesh point £ the diffusion-related terms in Eqs. (1.2) and (1.5) are modeled 

using interpolations that may involve the numerical values of the mesh points £ fl 2 (^i)- 
This contrasts sharply with the modeling of the convection-related terms which uses no 
interpolation. 


(b) In the a(2) scheme, the two sets of numerical variables associated with the individual 
meshes and fi 2 of the dual mesh are completely decoupled from each other. Contraxily, 
they are glued together in the a-p(I\) scheme through the interpolations referred to in (a). 


The conditions of flux conservation, Eq. (2.6) lead to the following equations, after 
integration and some algebra. Flux conservation for the element CE + (j, n), where the latter 
is defined, yields 


+ (“*)" 

-(a u]-( u,)” 

{* 


Ax — 

At 
2 


u n_1 
u j + 1 


/ \n— 1 Ax . 

- M j+ i — J Ax 


- n 


+ 

0. 


u n ~ l 
u j + 1 


/ \n-l At 
+ ( u t)j+l 2 


w(u x )] + tu(u*)" *] | A t 
w'{u x ) n ]+ 1 + } 


- h 


At 


(3.8) 


Flux conservation for the element CE_(j, n), where the latter is defined, yields 

,n— 1 


U 


” - K)” ^1 Ax - Ax 


+ {a - (u t ); y] - „ [w(u x r s + iK)-- 1 ] } At 
- (a u]:} + (u,)”:, 1 y] ~ i 1 + t»'K)?r 1 1 ]) At 


= 0. 


(3.9) 


It should be noted that during integration over the surface of CE + (j, n ), only h*(x, t ; j, n) and 
h*(x,t]j + l,n— 1) are used, while for CE_(j,n), only h*(x,t;j,n) and h*(x, t ; j — 1, n — 1) 
are used. Another point of interest is that the flux integration is made simple by the fact 
that the integral of a linearly varying flux over a straight line face of a CE is equal to the 
product of the length of the face and the value of the normal component of the flux at the 
midpoint of the face. 

In the interests of compactness and clarity, symbols are now defined for groups of quan- 
tities occurring frequently in the discretized conservation equations. Let 


v — a 


At 

Ax 




At 

Ax 2 ' 


(3.10) 
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(«:)” = T M" ; {<)" = Y <“*>" 


(3.11) 


The magnitude of u will be recognized as the Courant number, while £ is a mesh-related 
diffusion parameter that arises in computations involving diffusion equations. Eq. (3.4) can 
be recast as 

K)j = K)” + ! [K)” +1 - «);_,] (3-i2) 

for J > j > 0 and n > 0. Also, Eqs. (3.8) and (3.9) can be recast as 


-v L” - (<)” + £ [w(u+)] + w«)”~ 


+» «“+i + (««')J +1 1 - £ K«)"+i + «/«)£ 

= 0, J > j > 0, n > 0, 


(3.13) 


U 1 - (“!)“] - ["1-1 1 + 

+l/u"-(<)“ - f )” + w(u + )"' 


-V “"-I 1 + + w '( u t 

0, J > j > 0, n > 0. 


= 0, 


Collecting terms in Eqs. (3.13) and (3.14), we obtain 


(1 - v) u] + (1 + w£) (t*+) J - w'£ (t£)" +i + v (u + )” 

= (i - v) u*zi - (i - (<);;; - w£ («+)” 1 - v (ut )” 


(3.14) 


(3.15) 


(1 + u)u]-( 1 + w£) (t£)J + w'£ (<)”_ x - V (u+)” 

= U + ”) u l-i +{ l ~ (<)” 1 + ^ • 


(3.16) 


We turn next to a consideration of the equations at the initial line and at the boundaries. 
Based upon the specification of the initial values, Eq. (3.1), we set 


u° = Uj(xj), 0 <j < J. 


(3.17) 


The initial-value information is also utilized to specify the spatial derivative on the initial 
line; 

0<j<J. (3.18) 
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Similarly, based upon the specification of the boundary values, Eqs. (3.2) and (3.3), we set 

«o=ui(* B ). ”<>0, (3-19) 

u n J = u R (t n ), n > 0. (3.20) 

The boundary information is also utilized to specify the temporal derivative on the boundary: 

A t duj 


«): = 

= 


2 dt 
At du R 

2 dt 


(n, 


n > 0, 
n > 0. 


(3.21) 

(3.22) 


Eq. (3.12) and Eqs. (3.15)-(3.22) constitute three discretized equations per mesh point. 
On any fine of constant n, n > 0, one equation at each mesh point is given by Eq. (3.12) for 

0 < j < J, by Eq. (3.21) for 7 = 0 and by Eq. (3.22) for j = J. These equations can be used 
to eliminate (uf ) from the other two equations at each mesh point, as is done shortly. On 
the initial fine n — 0, the two remaining discretized equations per mesh point are given by 
Eqs. (3.17) and (3.18). On any line of constant n, with n > 0, the two remaining equations 
associated with each mesh point are given by Eqs. (3.15) and (3.16) for 0 < j < J, by Eqs. 
(3.15) and (3.19) for j = 0, and by Eqs. (3.16) and (3.20) for j = J. 

We proceed to eliminate (uf 'j using one equation, from the other two equations at each 
mesh point, in order to express tfie latter equations in terms of only the marching variables 
u n and {uf) n . At the boundaries, we will also substitute the boundary values given by Eqs. 

(3.19) and (3.20). It may be noted that the equations used to eliminate (iXf j , viz., Eqs. 
(3.12), (3.21) and (3.22), are applicable for any n > 0. Eqs. (3.17)-(3.20) do not involve 
a temporal derivative, thus leaving only Eqs. (3.15) and (3.16) to be transformed. When 

1 < j < J — 2, we may substitute for and (W") from Eq. (3.12) in Eq. (3.15) to 


obtain 


(i _ „) U J - (o+)"_ 1 + (l - * 2 + w() (<)’ + (j - »') « «)” +1 


(1 - v) - (l - V 1 - «/f) (<) 


+ 


(i- 


U) 




n— 1 


n — 1 


1 < 7 < J - 2 . 


When 2 < j < J — 1, we may substitute for (u'f'j and (uf 
(3.16) to obtain 

4 « 


n — 1 


(3.23) 

from Eq. (3.12) in Eq. 


(1 + i /)ti" 


= (1 + x) up - f + (l - X 2 - w'f) 


' j+1 
> n— 1 


(l - x 2 + w() (uj)" + (j + v) ( (»?)”_, 


+ 


(i + w )(( u *)" 1 


2 < 7 < J - 1. 


(3.24) 
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When 2 < j < J — 2, Eqs. (3.23) and (3.24) will constitute the interior equations of the 
a-pL(Il) scheme. 

At the left boundary j = 0, Eq. (3.19) is one equation associated with the mesh point. 
Substitution may be made in Eq. (3.15) for u" from Eq. (3.19), for (u+)” from Eq. (3.21), 

and for (u f'j from Eq. (3.12). This gives 

(i - V ) u L (n + a + «*) «)" - *>'( +‘'fw (n 

= (i-") - j (<)"-' - (i - s - *() (<))") 

+ ( j - w) $ (<)”“' , j = 0, (3.25) 


as the second associated equation. 

At the adjoining mesh point j = 1, one of the associated equations is Eq. (3.23). Substi- 
tution may be made in Eq. (3.16) for (u ~ )" from Eq. (3.12), for Uj_l from Eq. (3.19), and 

for 1 from Eq. (3.21). This gives 


(1 + U) U] - ^ (<)” +1 - (l - ^ + <) (<)" + (^ + 

= (1 + u)u L {t n ~ l ) + (l - «7f) (<)J_ 1 1 (u+)J 


,n-l 


+V 


A t duL 

~Y~dT 


(<”-) ■ 


3 = 1> 


(3.26) 


as the second associated equation. 

At the other boundary mesh point j = J, one of the associated equations is Eq. (3.20). 
Substitution may be made in Eq. (3.16) for u 7 - from Eq. (3.20), for (u + )” from Eq. (3.22), 

and for (uf ) . ^ from Eq. (3.12). This gives 


(1 + i/) u*(0 - (1 + wO (<)” + u/£ 


At duR 

dt 


2 
n— 1 


<<”) 


= (i + + (i - - “'d 

+ (j +5i ) 5 3=1 


(3.27) 


as the second associated equation. 

Finally, at the adjoining mesh point j = J - 1, one of the associated equations is Eq. 
(3.24). Substitution may be made in Eq. (3.15) for (u+)” from Eq. (3.12), for from 

Eq. (3.20), and for (ut'j from Eq. (3.22). This gives 


(i - v ) u i - t K)J_i + ( x - 1/2 + <) K)J + (i - w ') s K) 


j+l 
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(3.28) 


= (l-i/) u R (t n ! ) - (l - w'Z) (<)” +1 X “ < (<)" 


, n-l 


-v- 


A t du R 
2 dt 




j = J- 1, 


as the second associated equation. . 

It should be noted that Eqs. (3.18), (3.21) and (3.22) assume differentiability of the 
functions tt/, u L and u R , respectively. This differ entiability holds true for the numerical 
examples tested in this paper. However, a more general discretization of the initial and 
boundary conditions would be such that the discretized flux through the initial-line faces 
and the boundary faces of the CEs bordering the initial line and the boundaries equals the 
flux of the exact initial and boundary values, in an integral sense. For example, we would 


set 


and 


v2 - (u+)° Ax = [ 3 uj(x) dx, 1 < j < J, 

3 V 'i\ Jxj- Ax 


u j + {yZ ) j Ax = J ui{x ) dx, 0 < j < J - 1. 


(3.29) 


(3.30) 


For 1 < j < J- 1, Eqs. (3.29) and (3.30), rather than Eqs. (3.17) and (3.18), would together 
determine and (u+)°. The values of v? 0 and u°j are obtained from Eq. (3.17), while (uj) 0 
is obtained from Eq. (3.30) and (u+)° is obtained from Eq. (3.29). For the left boundary, 
we would set tn 

<-(u+)1 Ai= / u L (t)dt : 

L v /OJ Jt n -Ai 


n > 1, 


(3.31) 


and 


r / \ m rt n + At 

Uq + (u t + ) o At = J ui(t) dt , n > 0. (3.32) 

For n > 1, Eqs. (3.31) and (3.32), rather than Eqs. (3.19) and (3.21), would together 
determine Uq and (uT ] . The value of ^ is obtained from Eq. (3.32). The diffusive flux 
has not been accounted for in Eqs. (3.31) and (3.32), since it is not known from the exact 
boundary values. The right boundary values would be discretized in an analogous fashion. 
Such a general treatment would be valid even in the case in which the functions uj, ul and 
u R may be discontinuous. This completes the discussion of the treatment of the boundary 
conditions. 

Lastly, we express the interior equations, Eqs. (3.23) and (3.24), in matrix form with the 
aid of the following definitions 


T (j, n) 


def 





nW d ±f ( 0 -I/ £/ 4 \ 

V - 1 " ^ o 074 + <K ) 


QP * 


1 — 1 / 1 — I / 2 + w£ \ 

1 + u — (1 — i/ 2 + u>£) J 


(3.33) 

(3.34) 

(3.35) 
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n (l) def ( 0 (l//4 -W')£\ 

Qi - \ o -^/4 ; 


gC 2 ) dg 


0 0 

0 -i/f/4 


<?g ^ 


0 0 _ 
l + i/ 1 — i/ 2 — ty'£ 


Q o 

g(2) */ 


1 ( 2 ) dg ( 0 {v j 4 - ty)£ 

l 0 (i//4 + Ty)£ 


1 - v 

0 


(l — i / 2 — ty'£) 


0 


(3.36) 

(3.37) 

(3.38) 

(3.39) 

(3.40) 


Of = ( “ ^ /4 ) (3.41) 

With the preceding definitions, the interior Eqs. (3.23) and (3.24) may be expressed as 

Q-l~q(j - l,n) + Qo ] -q(j, n) + Qi ] ~q{j + 1, n) 

= Q-l~q (j - 2, n - 1) + Q^l ~q(j-l,n-l) + Q { q ] ~q(j,n- 1) 

+Q? V (j + l,n - 1) + Q 2 ] ~q(j + 2, n - 1), 


i.e., 2 

E Q, (1) T(j + ^)= E e! 2) T0 +/^-i). (3-42) 

Z=-l /=— 2 

3.2 The a-n(I2) Scheme 

This sub-section describes another implicit scheme, referred to as the a-ji(I2) scheme. In its 
specification, it differs from the a-/x(/l) scheme only in one respect, viz., Eq. (3.4) is replaced 

by 

(«.)? = -«(«.)? + + +1 - 2»D. J>}> 1- (3-43) 

This involves an alternative representation of the diffusion term in the governing equation. 
The flux conservation conditions for CE+ and CE_ are obtained as in the a-/z(/l) scheme, 
so that Eqs. (3.15) and (3.16) are also the other two equations associated with each interior 
mesh point in the a-/x(72) scheme. However, the expression for the time derivative occurring 
in these equations is that arising from Eq. (3.43) rather than that from Eq. (3.4). It may 
be noted at this point that in two special cases, the a-/x(il) and a-/z(/2) schemes become 
identical. The first case occurs when /x = 0, when the viscous terms disappear, and the two 
im plicit, schemes become identical to each other and to the a scheme described in Section 2. 
The second case occurs when a = 0. In the latter case, u = 0, and the time derivative terms 
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occurring in Eqs. (3.15) and (3.16) vanish because of the factor v. The two implicit schemes 
are then identical. 

Eq. (3.43) can be recast as 

j>]>\ (3.44) 


On substituting for ( ul ^ from Eq. (3.44) in Eqs. (3.15) and (3.16), we obtain the interior 
equations of the a-fi(I2) scheme: 


if “j-i + ( X “ v ~ if) u i + if u 

+ (l - v 2 + wZ) (u+)” - w'Z (ut). 

, (i 


+ f ) U ^-T U : 


3 + 1 

4 

- (i - ^ - s?«) (< )!’ - 


(3.45) 


and 


K„.n 






- (i - 1/ 2 + <) (u+)j + «/e (<)”_! 

7 U ; : 2 + ( 1 + z/_ y) u f Z l + 

(i - 1/ 2 - (y)"_! + (<)" 


n— 1 


. n— 1 


+ 


(3.46) 


The boundary conditions are discretized in the same way as for the a-/i(/l) scheme, i.e., 
Eqs. (3.17) through (3.22) are assumed. The conservation equations at or adjoining the 
boundaries may by obtained from Eqs. (3.15) and (3.16) by substituting for (t if) from Eq. 
(3.44) or from the known boundary values, in a manner similar to that used in the a-jj,(1 1) 
scheme. 

The interior equations (3.45) and (3.46) of the a-/x(/2) scheme may be expressed as Eq. 
(3.42), with the use of the following definitions: 


q(1) d g 


K/4 0 \ 

-K/4 w'Z ) 


n (l) def ( 1 -U- vZ /2 1 - V 2 + WZ \ 

Qo ~\l+u + uZ/2 -(1 -u 2 + wZ) ) 
n { 1 ) def ( VZI 4 -W'Z \ 

Ql - v -^/ 4 0 / 


(3.47) 

(3.48) 

(3.49) 
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q[ % t / 


( 0 °) 
{ vu 4 o ; 


(3.50) 


n (2) de/ / 0 0 _ 

^ 1 1 + v — i/£/2 1 — v 2 — u/£ 


(3.51) 


Qf H/ 


-^/4 \ 

i/^/4 y 


„(2) dej ( 1 - 1/ + ^/2 - (l - I/ 2 - w'£) 

V o o 


(3.52) 

(3.53) 


Qf d M 


-v £/ 4 0 
0 0 


(3.54) 


4 Solution Procedure 


We will first discuss the solution procedure for the a-fj,(1 1) scheme. The 2 (J + 1) equations 
represented by Eqs. (3.19), (3.20) and (3.23)-(3.28) relate the 2(J 4- 1) marching variables 
Uj and (u+)" at the nth time level to the 2 (J + 1) marching variables u” -1 and (ti+)” _1 at 
the (n — l)th time level. The equations thus define a time-marching scheme whereby the 
mar ching variables at the nth time level can be computed if the marching variables at the 
(n — l)th time level are known. Given values of the marching variables on the initial line, 
by Eqs. (3.17) and (3.18) for example, the scheme can be used to successively compute the 
marching variables atn = 1,2, 3,.... Next, we describe the solution of the implicit system of 
equations at any time level n. 

The flux conservation relations associated with CE + (j, n) , namely, Eqs. (3.25), (3.23) 
and (3.28), can be rewritten as 

(Eqs. (3.23) and (3.28)) 


and (Eq. 


(! - "K - f + 0 - 1,2 + w () (<)" + (l - ”') « (<)" +1 


= (S+)”- 1 , 1 < 3 < j ~ 1 

(3.25)) 

(1 + w?) (’4)" +l 

= j = 0 ' 


where 

(5+r 1 


— v) u’j+l 


4 




0 < j < J - 2 


(4.1) 


(4.2) 


(4.3) 


18 



and 


(s + )”-' = (i - v) ««(<”-■) - (1 


n— 1 
3 + 1 


j = J~ I- 


(4.4) 


Similarly, the flux conservation relations associated with CE_(j, n), namely, Eqs. (3.26), 
(3.24) and (3.27), can be rewritten as 
(Eqs. (3.26) and (3.24)) 

a + v) u] - £ (<)“ +1 - (i - s + w() (<)" + (j + •»') e 
= (S.)"- 1 , 1 < j < j - 1 

and (Eq. (3.27)) 

- (1 + w£) (u+)J + u/f 


(4.5) 


where 


+ 


+®) £ («£) 


\ n— 1 


71-1 

j 


and 


+w£ (u+)J 1 + V 


At duL 
2 dt 



j = j, 

(4.6) 

(1 -n 2 - 



'■<j<J 


(4.7) 

- ( u x, 

r 

>3 - 1 


- . 

3 = i- 

(4.8) 

are known from either the n — 

1 time 


level or from the boundary conditions have been transferred to the right-hand side. 

When 1 < j < J - 1, both Eq. (4.1) and (4.5) are applicable. Linear combinations of 
these equations imply that 


- ^ + (1 -v)w’ ?(<)”_, + 2(1-1 s + wt) (v,;)" + j-(i + v)w ' ; (uZ )“ 
= a + ^)(s + )r‘- (i-^xs- r 1 


j+1 

(4.9) 


and 


2 U " - w't ( u , + ) j +i - (<)”_,] = ( s *),- 1 + ( s -> r ' 


(4.10) 


It may be noted that the it" at the nth time level do not appear in Eq. (4.9). Note also 
that Eq. (4.10), obtained by adding corresponding sides of Eqs. (4.1) and (4.5), represents 
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flux conservation for the union of CE+(j, n) and CE_(j, n). Also, Eq. (4.10) implies that if 
w' = 0, then 

«“ = \ [(5 + )"-' + (S-)” -1 ] ( 411 ) 

The equations (4.2), which is applicable at j = 0, (4.9), applicable for 1 < j < J — 1, and 
(4.6), applicable at j = J, together form a complete system of linear equations for the (u+)” 
from which the u" at the nth time level are absent or are known from boundary values. After 
the (u+)” are computed by solving this system, the u 1 } in the interior can be found from Eq. 
(4.10). The u™ at the boundaries are known from Eqs. (3.19) and (3.20). 

The system of equations for the unknowns {u^) n - , j = 0, . . . , J is cast in matrix form 
with the aid of the following definitions : 

b 0 d = 1 + , c 0 d = -u/f , 

d 0 v ts + )r‘ - (i - *) ul (n - v Y^it (*”> ( 4l2 > 

a, = - [f + (1 - v)v/\ f 
bj 2(l-u 2 + wO 

c i d ~ [l “ C 1 + u ) w '} Z 

dj = f (1 + V) {S + T~ l -il-v) ( 5_)” _1 

aj d = w'£ , bj d = -(1 + wf) , 

^ *£ (s.)y 1 -(i + v )u R (n+-'Yw < - f "> (4 - 14) 


j = (4.13) 


v j ^ ( u x)j > j = 0,...,J (4.15) 

The system may be written in matrix form as Av = d, where v is a vector of length J + 1 
with ith entry v t , d is a vector of length J + 1 with ith entry dj, and A is a ( J + 1) x ( J + 1) 
tridiagonal matrix with subdiagonal, leading diagonal and supradiagonal formed with the 
entries a*, b t and c t respectively. Thus 



if j = i - 1 
if j = i 
if j = i + 1 
otherwise 


(4.16) 


It should be noted that the row and column indices of A and the row indices of v and d run 
from 0 to J. 

We next show that the matrix A is diagonally dominant in rows and columns when certain 
conditions are met. Let 

l>i/ 2 , £>0 (4.17) 
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and 


0.5 < w < 1 , w' = 1 — w 

Then Eqs. (4.12)-(4.14), the definitions of the entries of matrix A, imply that 

M = \bj\ = 1 + ^4 , |co| = M = (1 “ w )t 


\ a 3\ = 


(l -to) + -(2w- 1) 


\bj\ = 2(l -i/ 2 + w^j 


1/ 


(!-«») + y(2^-l) 


(1 - to) - -(2to - 1) 


M = 

where j = 1, . . . , J — 1. Thus 

N-K- 1 -tel > 2 

= 2 




1/ 


(1 — to) H- — (2w — 1) 


(l - v 2 + tu£) - (1 - to)^ - (2w - 1)4 


1 - H- (2to - 1) 


> 0 


where j = 1, . . . , J — 1- Obviously, we also have 

|6j| “ l a j+i I ~ l c i-i| > 0 i J = 2, . . . , J — 2 

Also 

N-N = N-|cj| = l + (2w-l)4>0 


|6o| — M > 1 + u>£ — 

= 1 + (2to — 1) 


\v\ 


(l-w) + -(2w-l) 




\bj\-\cj-x\ > 1 + — 

= 1 + (2 w - 1) 


4>o 


M 


(1-H + y( 2 ^-l) 




4>o 


|6i| — |co| — |o 2 | > 2(l-i/ 2 + <) -(1 -w)4- 


K 


(1 - to) + ■^•(2to - 1) 


= 2 


M' 


1 -i/ 2 + ( 2to- 1) | 1 -^U 


> 0 


(4.18) 

(4.19) 

(4.20) 

(4.21) 

(4.22) 

(4.23) 

(4.24) 

(4.25) 

(4.26) 

(4.27) 

£ 

(4.28) 
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I&j-il - |cj-a| - |«j| > 2 (l - v 2 + tof) - (1 - w)f - (1 - w) + y (2w - 1) £ 

= 2 1 - i/ 2 + (2tu - 1) ^1 - ^ f 

The above equations (4.19)-(4.29) imply that the ( J + 1) x ( J + 1) tridiagonal matrix A 
is strictly diagonally dominant in rows and columns. Thus the Thomas algorithm for solving 
the system is stable under the assumptions. 

If Eq. (4.17) is assumed to be true, and the assumptions of Eq. (4.18) are replaced by 

w > 1 and w 1 = 0 (4.30) 

then the entries of the matrix A reduce to 

N = M = 1 + , M = |co| = 0 



\bj\ = 2 (l — i/ 2 + u>£) , |Oj| = \cj\ = , j = -1 

From this, it can be shown that the matrix A is strictly diagonally dominant under the 
assumptions (4.30), and that the Thomas algorithm can therefore be used to perform the 
LU decomposition in a stable fashion. 

The LU decomposition of the matrix A and solution of the system is performed via the 
Thomas algorithm as follows. Let 

0o = (4-31) 

on = ■£- , 0i = bi- OLiCi-i , i = (4.32) 

Pi - 1 

This Gaussian elimination performs the LU decomposition of the matrix A = LU. The 
(J+l)x(J+l) bidiagonal matrix L has unity entries on the leading diagonal, and entries 
cti (i = 1, . . . , J) on the subdiagonal, where i is the row index, which goes from 0 to J. The 
( J + 1) x ( J + 1) bidiagonal matrix U has entries $ (i = 0, . . . , J) on the leading diagonal, 
and entries c* (i = 0, . . . , J — 1) on the supradiagonal. The system Av = LUv = d is 
transformed to Uv = L -1 d = <5 by a similar forward elimination process. The entries of 6 
are 

6 0 = d Q (4.33) 

Si = di — otidi-i , i = 1, . . . , J (4.34) 

The system Uv = 6 is then solved by a process of back-substitution 



(4.35) 
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and 


(4.36) 


Vi _ ( Q^z+l) i ^ 1 j • * * i 0 

Pi 

Finally, we mention that in the case of the a-fi(I2) scheme, prior elimination of the ti" at 
each mesh point to obtain a system for only the («*)" is not possible, as the conservation 
equations are implicit in both tt" and (u x )J. Hence, the resulting system for the u™ and the 
( w +) n must be solved by a block version of the Thomas algorithm described above. This 
makes the a-fx(Il) scheme computationally more efficient than the a-//(/2) scheme. 


5 Stability Analysis 

A rigorous discrete von Neumann stability analysis of the explicit ID CE/SE schemes for 
scalar conservation laws was carried out in [1]. In this section, we use the stability analysis 
procedure to analyse the a-n(Il) and a-/r(/2) schemes. For the rigorous development of the 
procedure, see [1]. 

We replace the discrete unknown with its discrete Fourier representation in the interior 
equation of the a-/i(1 1) scheme, Eq. (3.42). Then we use linearity of the interior equation 
to isolate a single Fourier mode. This is equivalent to making the substitution 

q{j,n) = <f(n,0) e^ e (i = V^T, -tt < 6 < tt) (5.1) 

in Eq. (3.42), where 9 is the phase angle variation per unit Ax of a single Fourier mode. 
Performing this substitution, we obtain 

<2 (1) (i/,£,ic,it/,6») <f (n, 6) = Q {2) (v,£,w,w f ,9) q*(n- 1,6), 

where 

Q m (v,£,w,w',6) 

= t^'el 11 
{ = - 1 

' i - i/ l _ „2 + W £ - if e ~ w + £ (* - w') e ie 
1 + v -(l-v 2 + wO-fe ie + t(z-w')e-' e 

and 


(5.2) 


(5.3) 


Q m ( v,S,w,ri,8) 

= 

l =- 2 

' (l-v)e' e -(l-^-a-wOO^-fe^ + C^ + w-l) 
(1 + v)e-' 6 (l-u 2 -(l -w’) 0 f + £ (f + 1 - w) 


(5.4) 
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We examine the determinant of to see when is invertible. We have 

A m (v,£,w,w',9) 

= det |Q (1) (u,£,w, w', 0)] 

= — 2 ^1 — i/ 2 ) + £(w — w'cos8) +i ^ — w'^j i/£sin0 . (5.5) 

From the above equation, we note that A^ ^ 0 for — 7r < 8 < n if 

u 2 < 1, £ > 0, and w > w 1 . (5.6) 

We assume that A (1) ^ 0. Then 1 exists, and we can multiply Eq. (5.2) from the left 

by [Q^J , to define the amplification matrix 

G{v,^w,w',8) = \Q {1 ) (v,£,w,w',8)} 1 Q {2) {v,£,w,w’,6). (5.7) 

An eigenvalue \(v, £, w, w 1 , 8) of G(v, £, to, w' , 8) is defined by 

G4> = \4> 

where </> is the corresponding (non-zero) eigenvector of G. Multiplying from the left by Q^\ 
and using Eq. (5.7), we have 

Q (2) <£ = AQ (1) 0. (5.8) 

Thus, any eigenvalue X(u, £, w, w ', 8) of G( v, £, w, to', 8) satisfies the condition 

det [Q (2) (i/,f,io,u/,0) - XQ w {v,£,w,w',0)] =0. (5.9) 

Using the definitions Eqs. (5.3) and (5.4), Eq. (5.9) is equivalent to 

AX 2 + B\ + C = 0, (5.10) 

where ^ 

A = (l — i/ 2 ) (to — to' cos 8) + i f - — to' j i/£sin0, (5.11) 


B = £ 1 — cos 8 + (to' — to) ( 1 + cos 8) — v 2 sin 2 
+iv 2^1 — v 2 ^ + (to + to' — 1) sin 8 


C = - (l - i/ 2 ) + £ [1 - to' - (1 - to)cos0] + i Q - to^ i^sin#. 


Thus, the eigenvalues of amplification matrix are given by 



(5.12) 

(5.13) 

(5.14) 
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These eigenvalues constitute the two amplification factors of the a-/r(/l) scheme. Note that 
there are two amplification factors rather than one, because there are two unknowns at each 
mesh point. Let 9 — 0. Then we find from Eqs. (5.20)-(5.14) that 


A 


+ 


1 an H \ = (1-i/ 2 ) 

i, and A_ ^ ^ + (1 _ ^ 


(when 6 = 0). 


Substituting 9 = 0 in Eqs. (5.3) and (5.4), and using the results in Eq. (5.8), we obtain 


0 —2 (1 — i' 2 ) 

0 2(1 -i/ 2 ) 


0 = 0 , 


whence we find that the eigenvector of G corresponding to A+ = 1 is given by 


$(9 = 0, A = 1) = a 


1 

0 


where a is an arbitrary constant. Based upon these facts, we term A + the principal am- 
plification factor, because it is exhibits the principal effect of the scheme on a wave of zero 
frequency. A_ is termed the spurious amplification factor. 

For general 9, the behavior of the amplification factors must be investigated numerically. 
Numerical experiments show that if either (i) 0.5 < w < 1 and w' = 1 — w, or (ii) w > 1 
and (re, w') lies between the lines w + w' = 1 and w — w' = 1 in the (re, w') plane, then the 
a-fi(Il) scheme is stable provided 

u 2 < 1 and £ > 0. (5.15) 


This condition on the stability applies, in particular, to the special case w = vJ = 1/2. 

Note that many other implicit solvers are unconditionally stable. However, the price paid 
for this “desirable” property usually is excessive numerical dissipation. Moreover, the use 
of a time-step size that is greater than that allowed by Eq. (5.15) generally results in a less 
accurate time-dependent solution. Thus we do not consider the more restrictive stability 
condition Eq. (5.15) to be a disadvantage of the a-fi{Il) scheme. 

We turn now to a closer examination of the amplification factors for the special case 
w = w' = 1/2. In Section 6, where the consistency and truncation error of the current 
implicit schemes are considered, it is shown that in this case, the a-/z(/l) scheme is second- 
order accurate in space and time. In this case, Eqs. (5.20)-(5.22) become 

J 4 = l-^ 2 + |(l-cos^) (w = w' = 1/2) 

B = £ [l — cos# — v 2 sin 2 + 2iv(\ - u 2 ) sin# {w = w = 1/2) 

C = — ( 1 — ^ 2 ) + I (1 — cos#) (w = w' = 1/2) 
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Then, corresponding to Eq. (5.10), we have 


-Q ± y/g 2 + (1 - n 2 ) 2 - £ 2 (1 - cosfl) 2 /4 
1 — v 2 + £ (1 — cos0) 


(5.16) 


where 

Q = - (l — cos 8 — v 2 sin 2 + iu(\ — v 2 ) sin#. 

In particular, when a = 0, i.e., when the scheme becomes a second-order accurate solver for 
the pure diffusion equation, then u = 0, and Eq. (5.16) can be simplified to yield 


1 — £sin 2 (#/2) 

+ = l + £sin 2 (#/2) ’ 


A_ = — 1 (w = w' = 1/2, v = 0) 


It is fo un d that A+ is identical to the amplification factor of the Crank- N icolson scheme ([28], 

p. 112). 

We next analyse the stability of the 2) scheme, with the same procedure as was 
used for the a-fi(Il) scheme. The substitution (5.1) in Eq. (3.42), with the definitions 
(3.47)-(3.54), lead to Eq. (5.2), where 

Q m (v,£,w,w' ,0) 

= E e'-'Q,' 1 ' 

i = - 1 

_ COS& + l - V - l - l/ 2 + w£ - w'£e t6 

— ^ cosd + 1 + v + *£ — (1 — v 2 + w£) + 

and 


Q {2) (v,£,w,w',8) 

= £ ellB Qi 2) 

i =- 2 

-^ + (1 - n + ^ )e ie - ^e 2ie - (1 - i/ 2 - (1 - w') f) e l6 - (1 - w) £ 
+ (1 + v - f)e~ ie + fe" 2 ^ (l-i/ 2 -(l-w')0<r* + (l-w)£ 

We examine the determinant of Q^\ to see when is invertible. We have 


(5.18) 


A (1) (i /,£,w,w',8) 


= det 

Q (1) ( v,£„w,w',8) 


= -2 

(l-n 2 )+e(a;- 

w' cos 8) — i ^1 + 8, sin 2 ^ w'v8, sin 8 


As in the case of the a-/x(/l) scheme, we find that A^ ^ 0 for — 7r < 8 < tt if Eq. (5.6) is 
satisfied. 
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We assume that ^ 0. Following the same logic as in the case of the a-/i(/l) scheme, 
we find that the eigenvalues of the amplification matrix for the a-/x(/2) scheme satisfy Eq. 
(5.9). This can be solved to yield the eigenvalues, as in Eq. (5.14), where Eqs. (5.20)-(5.22) 
are replaced by 

A= (l — i/ 2 ) + £ (w - w' cos 6) — i ^1 + £ sin 2 ^ w'v£ sin 6, (5.20) 

B = £ [1 — cos 6 + (w r — w) (1 + cos#)] 

T f 2 

+iv 2 (l - i/ 2 ) + (w + w' - 1)£ - (1 + w - w') (1 - cos#) sin0 (5.21) 

and 

C = - (l - v 2 ') +£ [1 - w’ - (1 - w) cos 9] +i ^1 - f sin 2 (1 - w)v£sin6. (5.22) 

Numerical evaluation of the amplification factors of the a-fi(1 2) scheme, shows that when 
(i) w > 0.5 and (ii) w' = w, the scheme is stable when Eq. (5.15) is satisfied. As pointed 
out in Section 3, when u = 0, the a-/i(/l) scheme and the a-/z(/2) scheme become identical. 
Hence, the a-p(/2) scheme shares with the a-/r(/l) scheme the property that when v = 0, 
its principal amplification factor becomes identical to the amplification factor of the Crank- 
Nicolson scheme. 

6 Consistency and Truncation Error 

In this section, we discuss the consistency and truncation error of, successively, the a-fi(Il) 
scheme and the a-ji{12 ) scheme. It will be indicated that, imder certain conditions, the 
interior discretized equations of the schemes are consistent with a pair of partial differential 
equations (PDEs). One of these PDEs is Eq. (1.2). In the following analysis, the grid 
point values of continuous functions of space and time are denoted by using grid indices as 
superscripts and subscripts. 

The flux conservation relation Eq. (3.24) is written for CE_(j + 1, n). This is symbolized 
by 

F-(j + l,n) = 0 (6.1) 

Eqs. (3.23) and (6.1) are two distinct conservation equations, although the conservation 
elements involved occupy the same physical region. The a-fi(Il) scheme can obviously 
alternatively be described by specifying these two equations. Linear combinations of these 
equations are next examined, in order to investigate the consistency of the scheme. Adding 
Eq. (3.23) and Eq. (6.1), and dividing throughout by 2 Ax At, results in 

[FDEl(u,u x )]”.~f =0 (6.2) 
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with FEE 1 as defined below. Also, subtracting Eq. (6.1) from Eq. (3.23), and dividing 
throughout by 2(1 — v 2 ) Ax 2 , results in 

[FDE2(u : u x )]"J =0 (6.3) 

Here, FDEl and FDE2 are defined as 

[FDEl(u, »)]“;| =' 

^ [«$» - «?«■ + «" - «r '] + [“"« - »“ + - «r ■] 


2 At L 

f!_L» 

2 Ax : ' +1 


; n 4 . ? , n - 1 _ ,, n - 1 

b + v j+i v j J 


p ( w ' + w — 1) ( Ax 
2 Ax h 4 At 


v j+i v ] 


v n ~ l -1- v n ~ l 
V j + 1 ^ V j 


L«-i 

+ 8Ax i J+a 


» - i;”- 1 + «; +1 - - v] +2 + vj - + vpi" 


[FPE2(n,t;)];;| ^ 


j [<?« + «" + + «rf 

+ 0/1^,^ [»*■» - »“ + »"+l - <-l + »“« - V T' + ““+1 1 - »#' 


+ 8 ( 1 ‘A, [«■+» - »" + »*+i - v u + «“« - -r 1 + "i+i 1 - v ti. 

~2 (i — t' 2 ) Ax [“"■« - ,j " + ““+> - u "~'] 

~ 2 (i-o ^ [“ ?+1 “ + " “ r ’] 

- 2 %Zo [*■ - « 5 Ti‘ + »" - »r‘] ( 6 - 5 > 

Let fi(x,t) and w(x,t) be smooth functions, and let 

p = v-— (6.6) 

Define, further, partial differential expressions PDE 1 and PDE2 by 

P W ^! + «g-,A (6.7) 

P B «(p)tf„ + ^| (6.8) 

1 n _ 1 

[FDFl(u, x)] . f and [FDF2(u, u)] f may be considered to be discrete approximations 

3+2 i J+ 2 

to [FDFl(u)] n 2 and [PDJ52(p)] n f , respectively. Then the errors ERl and ER2 in these 
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approximations are defined by 


[ Efli ]”;| =' [fdei(u, - ipdbhu))’‘;i ( 6 . 9 ) 

[BH2];;| =' [FDE2(u, S)]”;| - [P0£2(p)j;;| (6.10) 


With the aid of Taylor expansions, it may be shown that 


[ERl) n . + l = + — 


At -t~ 


d 3 u Ax 2 d 3 u At 2 


+ 


To~ 


dxdt ~ dx 2 dt 8 dt 3 24 
d 3 u Ax 2 d 3 u At 2 d 3 v Ax 2 d 3 v At 2 


+ CL 


dx dt 2 8 


8 dxdt 

It may similarly be shown that 

i ER Cl 

— v At 


dx 3 24 

+- S [a 2 At 2 — Ax 2 


- 




dxdt 2 8 


dx 3 24 

d 3 v At 2 / 3 \ 


du du d 2 u 
(1 — u 2 ) Ax | dt Q dx ^ dx 2 


£ (w — w ') 
z/£ Ax 


O (A 3 ) + 

d 3 v 


d 2 v Ax 2 
dx 2 8 


+ 


£ (w — w ') dv 
(1 — v 2 ) dt 
d 2 v At 2 
dt 2 8 


At 


24(1 -i/ 2 ) 

1 

24(1 - v 2 ) 
v At 


7 dx 3 ^ + 3 dxdt 2 At 

d 3 U 2 O a .2 

A * +3 &aP At 


dx 3 


dx 2 dt 


dt 3 


0( A 3 ) 


24 (1 — z/ 2 ) Ax 

■(3)0 M - (3)° ( A ‘) - (t3^° (**) 


(6.11) 


( 6 . 12 ) 


In Eqs. (6.11) and (6.12), all derivatives are evaluated at ((j + 1/2) Ax, (n — 1/2) At). Each 
symbol O (A 3 ) represents an infinite sum of terms, in which derivatives of u or v, and the 
quantities a, p and Ax l At m occur only as factors in the numerator of each term, with 
/, m > 0 and l + m> 3. 

Let u and v be a solution of the system of PDEs formed by PDEl{u) = 0 and p = 0. 
Then [F.D£l(fi,t5)]”~J and [FDE2(u, fi)]"~f are by definition the truncation errors of the 
discrete equations (6.2) and (6.3) satisfied by a solution of the a-p(Il) scheme. Also, now 
[PDEl(u)} n 3 ;l = 0 and [PD£2(p)]";| = 0. Hence, [FDEl{u, and 

[FDE2(u,v)] n ~l = [ER2} n ~f. Also, the first term on the right hand side of each of Eqs. 
(6.11) and (6.12) is zero. Assume i/ 2 ^ 1. Let the rule of mesh refinement be such that 
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A1 remains bounded as Ax — ► 0 and At —* 0. This implies that v and £ Ax also remain 
bounded. Examination of ERl and ER2 shows that the discrete equations (6.2) and (6.3) 
are in general consistent with only the stationary forms of PDEl(u) = 0 and PDE2(p) = 0. 
The second term on the right hand side of Eq. (6.12) represents a time-inconsistency, which 
makes the scheme only first-order accurate for u". If w = w', the a-p(Il) scheme is time- 
consistent, but is in general only first-order accurate in time because of the second term 
on the right hand side of Eq. (6.11). If, additionally, w + w' = 1, i.e. w = w' = 1/2, 
then the discrete equations of the a-p(Il) scheme are consistent with PDEl(u ) = 0 and 
PDE2{p ) = 0, and the scheme is second-order accurate in space and time. 

We turn next to a similar analysis of the a-p(I2 ) scheme. Manipulations similar to those 
leading to Eqs. (6.2) and (6.3) are performed on the interior equations (3.45) and (3.46) of 
the a-p(I 2) scheme. This yields the equations 

\FDEZ{u,u x )\ n -\ = 0 (6.13) 

[FDEA(u,u x )} n ~l = 0 (6.14) 

' 2 


Here FDE3 and FDEA are defined as 
[FDE3(u,v)) n } ;i = f 


1 


2 At 1 


u] +1 -u^l+u] 




vi+i-v'+vft-vr 1 


2 Ax L 

p (w ’ -f- w — 1) Ax 


2 Ax 


4 At 


ua 

T 


v"+i +«r 1 ' 


up 


■ 4 Ax2 [(«& - 2 up + «r‘) - - 2 up 1 + u]:l) - 

( u j+2 — 2u” + i — u ]‘j + ( u j+i — 2 uj + 


(6.15) 


[FDEi(u,v)\ n -\ tf 


»? + i + «; + <&? + <$ 


n— 1 




[(«;« - 2^ + -; + «r l ) + («& - 2»r’ + «#) + 


4 (1 — i/ 2 ) Ax 
(«”« - 2«? +1 + «?) + («?+, - 2< + «?-,)] 

1 


U j + 1 " + “j + 1 


_ u »-il 
i+i J 

<$» -«?«+<*" -«r‘] 


2(1 — i/ 2 ) Ax 
v 

2 (1 -i/ 2 ) Ax 

£ ( W - w/ ) U _ W n-X , n_ n-1 

2(1- „2) h+1 V i+1 +V i V J 


(6.16) 
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As before, let u(x,t) and v(x,t) be smooth functions, and let p be defined by Eq. (6.6)^. 
Then the errors ER3 and ERA in the approximation of [PDEl(u)] n .l and [PDE2(p)] ? 


by [FDE3(u,v)] n j+ ? and [FDE4(u,v)] n . + £ are defined by 


def 
i«~5 def 

b'+5 

With the aid of Taylor expansions, it may be shown that 
[ER3] n . 


j+i 


[Bfi3]” + | \FDE3(u,v))~i - [PDEUu)] 

[ERi ]";| =' {FDE4(u , «)]”;! - {PDE2(P)]"-? 


in-; 

■3+1 


dp . . d 2 v d 3 u Ax 2 


d 3 u At 2 




d 3 u Ax 2 


dx 3 24 
1 d 2 v 


+ d 


d 3 u At 2 cfti) Ax 2 


dx dt 2 8 


- A 1- 




dt 3 24 
cPv At 2 

dx dt 2 8 


8 dx dt 

It may similarly be shown that 

[ £ *Cf 

—v At 


a 2 At 2 - Ax 2 


dx 3 24 

d 4 u At 2 / 3 \ 

~ afJ, dx^di 4 +0 ( A ) 


(1 — z/ 2 ) Ax 
£ (w — w ') 

(1 - l/ 2 ) 

z/£ Ax 


<9tz du d 2 u 

dt a dx ^ dx 2 

d 2 v Ax 2 


C>( A 3 ) 

d A u 


+ dx 2 8 + 8 


£ (io — «/) 
(l-i/ 2 ) dt 
d 2 t) At 2 


At 


24(1 - I/ 2 ) 

1 

’24(l-i/ 2 ) 
z/ At 


ax 4 

a 3 u 


A 2 O A j.2 

Ax + 3 „ ^ At 


dx 2 dt 2 
d 3 u 


dx 3 ^ + 3 dx dt 2 


At 2 


n d 3 u 2 d 3 u 2 

3 d^di Ax + dP A 


24 (1 — I/ 2 ) Ax 
z^A^ / A s\ I o f A 3\ yA f_o (A 3 ) 

+ (l_ I/ 2) U l A j (! — Z/ 2 J V ) (l-z/2)Ax V / 


+ o (a 3 ) 

i/ At 


(6.17) 

(6.18) 


(6.19) 


( 6 . 20 ) 


In Eqs. (6.19) and (6.20), all derivatives are evaluated at (( j + 1/2) Ax, (n — 1/2) At). 

It is seen that the expressions for ER3 and ERA axe very similar to those for ERA and 
ER2. The same remarks made regarding the consistency and truncation error of the a-p(Il) 
scheme are applicable to the a-p(I2) scheme. 

The main conclusion of this section is that if (i) v 2 ± 1, (ii) the rule of mesh refinement 
is such that ^ remains bounded as Ax — * 0 and At — > 0, and (iii) w — w 1 — 1/2, then 
the interior discrete equations of the a-p(Il) and a-p(I2) schemes are consistent with the 
advection-diffusion equation PDEl(u) = 0 and with PDE2(p) = 0, and the schemes are 
second-order accurate in space and time. 
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7 Filtering of the Initial Conditions 


The accuracy of the solution may be improved by a modification of given initial conditions. 
T his modification represents a filtering out of spurious components from the numerical rep- 
resentation of given initial conditions. 

Let A t = 0, i.e., u = £ = 0. Then Eqs. (4.9) and (4.10) of the a-n(Il) scheme become, 
for j = 1,..., J - 1, 


(o;=i[us; -(C)] 


and 


“’ = 5 




Eqs. (7.1) and (7.2) can be rewritten as 


. n—1 


5 [« + « 


and 


l 

4 L 


+ u 'i-i + 2 < 1 - (“j)” +1 ‘ + 


(7.1) 

(7.2) 


(7.3) 


(7.4) 


where the terms on the right hand sides of Eqs. (7.1)-(7.4) are all known from the initial 
time level, given by Eqs. (3.17) and (3.18). 

Corresponding to Eqs. (7.3) and (7.4), the conservation equations (3.15) and (3.16) at 
the boundaries can be written as 


and 


1 

2 

1 

2 


[«+«r 


1 

2 


- «" + 


- «£? + 


(“C - 

K 

v n— 1~ 

'j+1 J 

, 3=0 

(7.5) 

«r - 

« 

ri 

'j-ij 

, 3 = J 

(7.6) 


where the right hand sides of Eqs. (7.5) and (7.6) are known from the initial time level or 
the boundary values. 

Eqs. (7.1), (7.2), (3.19), (3.20), (7.5) and (7.6) serve to determine ‘filtered’ initial val- 
ues. The von Neumann stability analysis of the interior filtering ‘scheme’ is performed as 
previously, with the substitution of Eq. (5.1) to obtain the form of Eq. (5.2) where 


Q w = 


1 1 
1 -1 


and 


Q< 2) = 



(7.7) 
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Thus [<2 (1) ] 1 = ^Q (l) and therefore the amplification matrix of the interior filtering scheme 
is given by 

G = [q (1) ] _1 Q (2) = 


cos 9 —i sin 9 
i sin 9 — cos 9 


(7.8) 


8 The Dual- Mesh Explicit <;-// Scheme 

In this section, we present an explicit numerical scheme to solve Eq. (1.2). It differs from 
the explicit a-fj. scheme presented in [1, 6] in that the diffusion term in Eq. (1.2) is better 
modeled in the current explicit scheme, and the diffusion flux in Eq. (3.6) is modeled more 
flexibly. These differences both introduce a coupling between the meshes fti and 0 2 . This is 
in contrast to the explicit a-p scheme of [1, 6], which is formulated on fli (or fl 2 ) exclusively. 
Hence, the present explicit scheme will be termed the dual-mesh exphcit a-p scheme. 

The exphcit scheme will be presented by indicating the modifications to the a-fi(Il) 
scheme to ch an ge it to an explicit scheme. The dual space-time mesh, the SEs, and the CEs 
are defined as for the a-fi(Il) scheme. The numerical representation defined in Eq. (2.1) 
is also assumed. We continue to take advantage of the notation defined in Eqs. (3.10) and 
(3.11). Instead of using Eq. (3.4), the diffusion term in that equation will be lagged. Thus, 
instead of using Eq. (3.12), the unknown («£ + )" at the time step n will be expressed as 



0 < j < J- (8.1) 


However, the known quantity ( j at the time level n — 1 continues to be calculated from 

Eq. (3.12), rather than from Eq. (8.1). Calculation of (<)” 1 from Eq. (8.1) would result 
in a three-time-level scheme, which would require additional computer memory and would 
require knowledge of the solution at n = 1 in addition to the initial condition at n — 0. 

The definition of the diffusion flux is also altered by setting w' — 0 in Eq. (3.7), so that 
for (x,t) E SE(j,n), 


ic / , . \ del 

U x {x,t]j,n) = 


w(u x )^ + w(u x )™ 1 , if t < £ n ; 
(u*)7 , if t > t n . 


( 8 . 2 ) 


where w is a weighting factor to be specified, and w d = 1 — w. This definition is used to 
complete the definition of the numerical space-time flux in Eq. (3.6). The numerical flux 
analogue is required to satisfy conservation over the CEs. Thus, Eqs. (2.6) are required to 
be satisfied. With the modifications described above, this leads to Eqs. (3.13) and (3.14) 
being replaced by 



v n ~ l 
u i + 1 



n— 1 

j + 1 
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-v u"-(u ( + )J +f [w«)” + u;«)" l ] 

+!/ U^ 1 + (u+ )” +i : ] - e(0”+l 

= 0 , 0 < j < J (8.3) 


and 



= 0 , 0 < j < J. (8.4) 

As with the implicit schemes, the discretized initial and boundary information required by 
the scheme will be assumed to be given by Eqs. (3.17)-(3.22). However, the same remarks 
as before apply more general treatment of the initial and boundary conditions. Collecting 
terms in Eqs. (8.3) and (8.4), we obtain 

(i - u) u? + (i+ (t£)j + v (u+y. 

= a—*') - (i - o - n K)” -1 - * K )"~! (8-5) 

and 

(1 + v) u? - (1 + w( ) (tt * )” - V (<)” 

= (1 + v) u]:i + (1 - {) («+)""* + W( (uj)” _1 + v (<)”;' . (8.6) 

Eq. (8.5) is valid for 0 < j < J . Together with the boundary condition Eq. (3.19), it 
supplies one equation per mesh point ( j,n ), 0 < j < J. Eq. (8.6) is valid for 0 < j < J 
at the time level n. Together with the boundary condition Eq. (3.20), it also supplies one 

equation per mesh point ( j,n ), 0 < j < J. The quantities (uf ) ^ and (uf appearing 

in Eqs. (8.5) and (8.6) are known from the previous time level n — 1, through Eq. (3.12), 
or from the boundary data Eqs. (3.21) and (3.22). The quantities (ut) j at the current 

time-level n will be eliminated from Eqs. (8.5) and (8.6) as follows. 

For j = 0, the known boundary data from Eqs. (3.19) and (3.21) are used in Eq. (8.5) to 
yield 

(1 + »?) («+)" = (5+)’- 1 . 7-0, (8.7) 

where 

(sur 1 = d - *> - d - o (“t)”"' - " K 

-a 
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Eq. (8.7) and the boundary condition Eq. (3.19) constitute two equations at j = 0. 
For 1 < j < J — 1, we substitute for (ut) . from Eq. (8.1) in Eq. (8.5) to obtain 


(1 - v)uj + (l - v 2 + wt) (<£)" = (5+)” 1 


( 8 . 8 ) 


where 

Similarly, for 1 < j < J — 1, we substitute for from Eq. (8.1) in Eq. (8.5) to obtain 

(l + n)«;-(l-n 2 + m{)(<)” = ( 5 -)r 1 < 8 ' 9 ) 


where 


(s-)"- 1 = (i + *o ”"- _ i 1 + (i - « - t ) «: + 




\ 71 — 1 


i+i 


+ 1/ 




Eqs. (8.8) and (8.9) constitute two equations at each mesh point (j, n) such that 1 < j < 
J —l. 

Finally, at j = J, the known boundary data from Eqs. (3.20) and (3.22) are used in Eq. 
(8.6) to yield 

- (i + tuO (t£)J = (s_)" -1 , j = J, ( 8 - 10 ) 


where 


(S-)- 


71 — 1 


(1 + u) U# + (1 - 0 (<)”_' + < (<)” 

- (i + v) u" + u (uty. 


+ V 



Eq. (8.10) and the boimdary condition Eq. (3.20) constitute two equations at a boundary 
mesh point with j = J. 

All mesh variables on the right-hand sides of Eqs. (8.7), (8.8), (8.9) and (8.10) are known 
from the marching variables at the previous time level or from boundary values. The system 
of equations described above, with two equations per mesh point, enable a marching step in 
which the solution vR and (u+)" at the current time level n can be determined by explicit 
solution. The explicit solution process will now be discussed. At the left boundary j = 0, u" 
is known from the boundary condition, Eq. (3.19), while )j is obtained from Eq. (8./) 
as 

(<)“ = (s + )"-‘ / (1 + o l 8 ' 11 ' 
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At the right boundary j = J, u" is known from the boundary condition, Eq. (3.20), while 
(u+)” is obtained from Eq. (8.10) as 

(t£)” = -(S_)”-7(l + u,fl (8.12) 


At a mesh point (j, n) with 1 < j < J — 1, u" and (u+)" are obtained by solving Eqs. (8.8) 
and (8.9) simultaneously to yield 


1 r. 


^ = 2L 


(s + )r + (s-)” 


n— 1 


(8.13) 


and 



1 

2(1 — v 2 + w£) 


(1 + v) (S + )”-‘ 


(i - V ) (5.);- 1 ] . 


(8.14) 


Thus the solution at the new time level can be explicitly calculated. Note that Eq. (8.13), 
which does not involve (u+)", is obtained by simply adding together the conservation equa- 
tions (8.8) and (8.9). It represents flux conservation for CE(j, n), the union of CE+(ji, n) and 

CE_0». 

The following important facts should be noted. The boundary information is used in the 
dual-mesh explicit scheme in the same way as it is used in the a-fi(Il) scheme, i.e., the 
boundary equations of the two schemes are the same. If the conditions 


«? = «r 1 




(8.15) 


hold for 0 < j < J and some n — N, then the remaining equations of the dual-mesh 
explicit a-fi scheme written at the time level N become identical with those of the a-fi(Il) 
scheme written at the same time level. Specifically, Eqs. (8.1) and (3.12) become identical, 
Eqs. (8.2) and (3.7) become identical, and therefore the integral conservation relations Eqs. 
(8.5) and (8.6) of the dual-mesh explicit scheme become identical with those of the a-/i(/l) 
scheme, Eqs. (3.15) and (3.16), respectively. Hence, solution time-slices u” and u” _1 that 
satisfy the condition (8.15) and the discretized equations of the dual-mesh explicit scheme, 
will also satisfy the discretized equations of the a-/x(/l) scheme. 

Next, the stability of the dual-mesh explicit scheme will be discussed. For 1 < j < J — 1, 
we substitute for from Eq. (8.1) and for (uf ) ^ and ^ from Eq. (3.12) in 

Eqs. (8.5) and (8.6) to obtain the interior equations of the scheme in terms of u" and (u^)", 
i.e., 


(l-i/)u”+(l -^ 2 + <) (t£)J 

0- - "Ks- - t ("C-R-^-s)] 


(8.16) 
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and 


(1 + v) u n j - (l - i/ 2 4- tof) (m+)” 

(1 + V ) f (<);;; + [i - - 2 - « (i + j)] «: 


+ 1 J+“ 


U 


T + tM 


n— 1 

i+i 


(8.17) 


Let 1 - i/ 2 + u>£ ^ 0. Then, as is shown by Eqs. (8.13) and (8.14), the interior equations 
(8.16) and (8.17) can be solved for tt" and (ti+)”. Let ~q (j, n) be given defined by Eq. (3.33). 
Then this solution is written in the form 


V 0, n ) 

- Q-^~q (j — 2,n — 1) + Q- 2 ~q(j - 2,n - 1) + QoVO> n - 1) 

Qi VO + l) n — 1) + Q 2 VO + 2,n — 1) 


i.e. 


where 


V 0 , n ) = J2 Qi~qU + l , n ~ 1 ) 

f =-2 




0 

0 


<fe / 1 

V-l — o 


1 + 1 / 


I 

4(1— i/ 2 +iu£) / 

1 — i/ 2 — £ 


2 ( i=£ \ c ( 1 ~*H 1 ~ l,) ( 1 - /a ) 

\ \1— i/ 2 +w£ ) 1— i/ 2 +iu£ 

( 0 




i £ 
2 


(8.18) 

(8.19) 

( 8 . 20 ) 

( 8 . 21 ) 
( 8 . 22 ) 

) (8-23) 

4(1— i/ 2 -fti;£) / 

The stability analysis procedure is the same as that used in Section 5 for the implicit 
schemes. When the substitution of Eq. (5.1) is made in Eq. (8.18), the result is 


0 


Qi = - 
V1 2 


2c(»-i+4) 

1— l/ 2 +U>£ 

-(i-i/ 2 -0 


1 — 1/ 

i-„» €(i+f)-0+^)(i -^) 

1 — l/ 2 4*u;£ 1— 




0 

0 


-k(i+i/) 


r(n,0) = Q(i/^,^,^)r(^-i^). 


(8.24) 


Here, 



/=— 2 


<2n Q 12 
Q21 Q22 


(8.25) 
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where 


(8.26) 

(8.27) 

(8.28) 


Qn = cos# — ivsmQ 
Q\2 — ~ sin 2 # — 2i (l — v 2 — sin#j 
Q 21 = i ^1 — v 2 ^ sin# / (l — v 2 + w^j 

022 = 2 (1 - ? T <) ( " cos2tf + isin2e) + * (“ ~ 1 + t) 

-2 (l - i/* - f) cos# + i [i/f - 2i/ (l - v 2 )] sin#j (8.29) 

Examination of Eq. (8.24) shows that Q( u, £, w, #) is the amplification matrix for the interior 
scheme. The eigenvalues of Q are the amplification factors. An eigenvalue A of Q satisfies 
Eq. (5.10) with 

A = 2(l-*/ 2 + <) (8.30) 

B — — £ 2(1 + w)cos# + ^ 2 sin 2 # + 2(w — l)j 

— i^sin# (1 — 2w)£ — 4 ^1 — i/ 2 ) — £ cos#j (8.31) 


C = 2(w — 1)£ cos # — sin 2 # — 2 (l — i/ 2 + w^j 

+ii/£sin0(l — 2u; — cos#) (8.32) 

The eigenvalues X± have been evaluated numerically for a range of values of v, £, w and #. 
These numerical experiments show that if 

w > 2 (8.33) 

then the dual-mesh explicit scheme is stable provided only that Eq. (5.15) holds, i.e., if 
v 2 < 1 and £ > 0. 

We note the unusual fact that if Eq. (8.33) is satisfied, the magnitude of the viscosity has 
no influence on the stability condition of the dual-mesh explicit scheme. A similar statement 

can be made about the explicit a-/i scheme of [6]. This is in contrast to most explicit 

schemes for diffusion equations, for which the value of £ is limited by stability constraints. 
This fact, together with the remarks made earlier concerning the condition (8.15), points to 
an important possibility for the dual-mesh explicit scheme. The dual-mesh explicit scheme 
ra.n be used to obtain ‘steady-state’ time-asymptotic limits of solutions to the convection- 
diffusion equation in an efficient manner. Time steps as large as those used with the a-fi(Il) 
scheme, and limited only by the CFL criterion, can be used. At each time level, fewer 
airthmetic operations are required for the explicit scheme than for the 1) scheme. This 
computational advantage would be even greater when multidimensional extensions of the 
algorithms are considered. The steady-state solution time-slice will be the same as that 
obtained with the a-ji(Il) scheme, and thus be second-order accurate in space. 


38 



9 Numerical Results 

Three test problems will be used to evaluate the accuracy of the implicit a-fi schemes. The 
weighting factors for the viscous fluxes are taken to be w = w' = 0.5 for all the computations 
run with the implicit schemes. Unless stated otherwise, the a-/x(/l) scheme was used for the 
computations in subsections 9.1 through 9.3. 

9.1 Decaying Traveling Sine Wave 

In the first problem, we consider a special case of Eq. (1.2) (a = 1 and p = 0.01) with 
0 < x < 1 and t > 0. The initial and boundary condition functions uj(x), UL{t) and u^{i) 
are defined such that they are consistent with a special solution to Eq. (1.2), i.e., 

u = u e (x, t) = f exp(— 47 r 2 fj,t) sin [27 t(x — at)] (9.1) 


At any time t = t n , let 


Li{u) = 


j- i 


(J - 1) exp(— 47r 2 /at n ) fr[ 


Y K-^.m 


(9.2) 


L\{u) is an error norm (per mesh point) which is normalized by the decay factor of u e (x, t). 

Let J — 80 (i.e., Ax = 1/80) and v = 0.8 (i.e., A t = 0.01). Then L\{u) = 0.1641 x 10 -3 
at t — 4 (i.e., n — 400). Through numerical experiments, it has been shown that Li(u) at a 
given time t is reduced by a factor of 4 if both Ax and At are reduced by half, confirming 
the second-order accurate nature of the scheme. 

Fig. 8 shows a comparison of the computed solution at t = 4 with the exact solution. It 
also shows the error u e (xj,t n ) — u™ scaled by the peak magnitude of the exact solution at 
that time level. It is seen that the maximum error is less than 0.3% of the peak magnitude. 
The peak magnitude is seen to be just over 0.2. 

Fig. 9 shows a comparison of the errors in the solutions obtained with the a-fx(Il) scheme 
and the implicit MacCormack scheme ([28]). Both schemes were applied with the same 
parameters and mesh as described above. It is seen that the current scheme is considerably 
more accurate than the implicit MacCormack scheme. Refinement of the grid keeping v 
constant would actually make the comparison even more favorable. This is because the 
implicit MacCormack scheme is second-order accurate in space and time only if f is held 
constant and v — * 0 when refining the grid. 

9.2 Pure Diffusion 

We consider a special case of the convection-diffusion equation with a = 0 and p = 1, in 
the domain 0 < x < 1 and t > 0. The initial/boundary conditions completing the problem 
specification are (i) u(0,t) = u(l,f) = 0 for t > 0, (ii) u(x, 0) = 2x for 0 < x < 0.5, and 
(iii) u ( x ,o) = 2(1 - x) for 0.5 < x < 1. The solution u(x,tj exhibits the diffusive decay 
of the initial sawtooth shape. An exact series solution is available, see for e.g. p.15 of [29]. 
For the CE/SE computation, uniform mesh intervals Ax = 0.02 and At = 0.005 are used. 
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Fig. 10 shows the time-slice at t = 0.05, comparing numerical and exact solutions, and also 
showing the error scaled with the peak exact value at that time level. The maximum error 
magnitude is seen to be about 0.5% of the peak solution value. At t = 1 (not shown), when 
the peak solution value has dwindled to about 4 x 10 -5 , the maximum error magnitude is 
about 0.15% of the peak solution value. 

The same problem on the same mesh was solved with the Crank-Nicolson scheme, which 
is probably the best traditional scheme for the parabolic pure diffusion equation. Fig. 11 
shows a comparison between the current scheme and the Crank-Nicolson scheme, at t = 0.05. 
The schemes are seen to be of similar accuracy, as is to be expected from the identity of the 
amplification factors. The current scheme is seen to be slightly more accurate for this case. 

9.3 Steady State Boundary Layer 

We next consider the problem defined for the convection-diffusion equation in the domain 
0 < x < 1 and t > 0 by the conditions (i) u(0,t) = 0 for t > 0, (ii) u(l,t) — 1 for t > 0, and 
(iii) u(x, 0) = x for 0 < x < 1. The ‘steady-state’ or time-asymptotic limit of the solution is 
u(x,oo) = [exp (ax///) — 1] / [exp(a//t) — 1]. 

The case a = 1, /t = 1 (i.e., ‘Reynolds’ number Re = a//x = 1) is first considered. Fig. 12 
shows the steady-state limit of the numerical solution, compared with the exact steady-state 
limit. With only three interior mesh points, the numerical solution has a maximum error of 
about 0.001. 

In Fig. 13, we compare the errors in the a-//(/l) and a-/t(J2) solutions for the same 
problem, but with \i = 0.1, so that Re = 10, which leads to a steady-state boundary layer 
at x = 1. Mesh spacing Ax = 0.05 was used. The two schemes yield very similar results. 

The final case considered is the same problem, but with fj, = 0.01, so that Re = 100. 
This leads to the formation of a fairly sharp boundary layer, because the thickness of the 
layer scales as the inverse of the Reynolds number. Uniform mesh intervals Ax = 0.0025 
and At = 0.002 are used, so that the Courant number is 0.8. Fig. 14 shows the computed 
and exact steady-state limits, together with the error. The boundary layer is seen to be well 
resolved, with the maximum magnitude of the error being about 1% of the solution peak. 

9.4 Calculations with the Dual-Mesh Explicit Scheme 

In this subsection, we examine results obtained with the dual-mesh explicit a-fi scheme for 
the same problems as in subsections 9.1 and 9.3. The value of the weighting factor w in the 
dual-mesh explicit scheme was always taken to be 2. 

Fig. 15 shows the solution time-slice at t = 4, computed with the dual-mesh explicit 
scheme, for the same problem as considered in subsection 9.1. The same spatial and temporal 
mesh intervals were used as in subsection 9.1, so that Fig. 15 can be compared with Fig. 8. 
It is seen that the a-// (71) scheme is considerably more accurate than the dual-mesh explicit 
scheme, owing to the latter being only first-order accurate in time. The peak error in Fig. 
15 is just over 7% of the solution peak, as compared to only about 0.2% in Fig. 8. The result 
obtained for the same problem with the same mesh parameters, using the explicit 1) 
scheme of [6] is shown in Fig. 16. The error plotted in Fig. 16 is very similar to that in Fig. 
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15. This shows that the error is dominated by that arising from the first-order-time-accurate 
modeling of the viscous fluxes in the integral conservation equations in both explicit schemes. 
Fig. 17 shows the solution time-slice at t = 4, computed with the dual-mesh explicit scheme, 
but using a time step that is four times smaller than that used for Fig. 15, while keeping the 
spatial mesh intervals the same. The error is smaller by about a factor of ten than that in 
Fig. 15. Numerical experiments show that the numerical space-time solution can be made 
to approach the exact solution as closely as desired by refining the mesh with a refinement 
rule that holds £ constant rather than v, so that At is reduced proportionally to Ax 2 rather 
than to Ax. The same is true of the explicit a-/x(/l) scheme of [6]. For multidimensional 
flows, the explicit schemes could have an advantage over the implicit a-fi(Il) and a-/i(/2) 
schemes, with regard to computational cost to achieve a given level of accuracy. 

The dual-mesh explicit scheme was also applied to the ‘boundary layer’ problem of sub- 
section 9.3. Fig. 18 shows the ‘steady-state’ limiting behavior of the numerical solution, at 
t = 3. The computation was performed with the same values of Ax, At, v, and £ as were 
used for the a-fi(Il) scheme in subsection 9.3. The same initial and boundary conditions 
were used. From Figs. 14 and 18, and from an examination of the computer printout of 
the values, the steady-state limit reached by the a-fi(Il) and dual-mesh explicit schemes are 
the same. This is to be expected, on the basis of the remarks at the end of Section 8. Fig. 
19 shows the ‘steady-state’ solution time-slice obtained with the explicit a -// scheme of [6], 
using the same conditions as for the other schemes. There is seen to be a considerable error 
in the steady-state limit, owing to the lack of modeling of physical diffusion in the discretized 
differential form of the conservation law. This error can be reduced to match the error level 
in Fig. 18, at the cost of increased computation, by reducing At and perhaps reducing Ax. 

For multiple spatial dimensions, the computational cost of an implicit time-step is gener- 
ally much higher than that of an explicit time-step. The analysis of the dual-mesh explicit 
scheme, confirmed by the numerical results, shows that the scheme can be used to obtain the 
same steady-state limit as with the a-fi(Il) scheme. It is expected that for most steady-state 
results, the steady state can be reached in approximately the same number of time steps by 
the two schemes. Hence, when generalized to multidimensional flows, use of the dual-mesh 
explicit scheme may be the more efficient way to obtain the steady state with a desired level 
of accuracy. 

10 Summary and Conclusions 

In the explicit a-p. scheme [6, 4], the diffusion term in Eq. (1.2) is not modeled, i.e., Eq. (2.3) 
is assumed. Also the diffusion term in Eq. (1.5) is modeled with no interpolation or extrap- 
olation, with a resulting reduction of time-accuracy. Obviously, such a solver can be used 
only when the diffusion term is small compared with the convection term. Contrarily, the 
diffusion terms in both Eqs. (1.2) and (1.5) are modeled in the current implicit solvers. The 
diff usion term in Eq. (1.2) is also modeled in the present dual-mesh explicit a-fi scheme, 
although it remains first order accurate in time because of the explicit nature of the modeling 
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of the diffusion fluxes in Eq. (1.5). 

The current implicit and explicit solvers are carefully constructed so that they become 
identical to the explicit a scheme for the pure convection equation, when the viscosity coef- 
ficient vanishes. Because the a scheme is nondissipative, this construction ensures that the 
physical dissipation is never overwhelmed by numerical dissipation in the present implicit 
and explicit solvers. 

The implicit schemes have been shown to be stable provided the Courant number does not 
exceed unity in magnitude, for certain ranges of the weighting parameters w and w' . In these 
cases, there is no dependence of the stability of the schemes on the viscosity parameter. This 
is true, in particular, for the case w = w' = 1/2. Similar remarks concerning the stability 
are true for the dual-mesh explicit a-\x scheme, for values w > 2. Stability analysis reveals 
the remarkable facts that (i) if fj. = 0, the amplification factors of the dual-mesh implicit 
and explicit a-ji schemes reduce to those of the classical Leapfrog scheme, and (ii) if a = 0, 
one of the amplification factors of the implicit a-p schemes reduces to that of the classical 
Crank-Nicolson scheme, while one of the amplification factors of the dual-mesh explicit a-/j 
scheme reduces to that of the DuFort-Frankel scheme. 

The truncation error analysis of the discretized equations of the implicit schemes shows 
that in general they are consistent with the convection-diffusion equation, and are first-order 
accurate in time. For the case w = w' — 1/2, which will be the case used in practice, the 
scheme is second-order accurate in time and space provided the Courant number remains 
bounded when re fining the space-time mesh. Although a truncation error analysis of the 
dual-mesh explicit a-ji scheme was not performed, it is evident that the scheme is second- 
order accurate in space and first-order accurate in time. 

Numerical examples have borne out the conclusions of the stability and truncation error 
analysis. The a-fi(Il) and a-pt(1 2) schemes have been shown to be similar to each other in 
their properties and performance. The current implicit schemes were seen from the numerical 
examples to enable stable accurate computations over the whole viscosity range, from the 
pure diffusion case to convection dominated problems. The dual-mesh explicit a-ji scheme 
was shown to be useful for the same range of problems, though needing a smaller time step 
to achieve comparable accuracy in the time-dependent problems. In the problems where 
only a steady-state limit is of interest, the dual-mesh explicit a-ji scheme yields the same 
spatial accuracy as the implicit a-p{I\) scheme, with potentially much lower computational 
cost when extended to problems in multiple spatial dimensions. 

The analysis and results noted in this study indicate that the development of the dual- 
mesh implicit and explicit solvers along similar fines for the time-dependent compressible 
Navier-Stokes equations may prove equally fruitful. 
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Figure 3 — The space-time mesh of the a scheme 
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Figure 6 — The space-time mesh of the implicit schemes 
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Figure 9 — Sine Wave Example : Comparison of Schemes 
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Figure 1 1 — Pure Diffusion : Comparison of Schemes 
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Figure 12 — a-\i (II) Scheme : Boundary Layer, Re = 1 



Figure 13 — Boundary Layer, Re = 10 : 
Comparison of Schemes 
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Figure 14 — a- |X (71) Scheme : Boundary Layer, Re = 1 
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Figure 16 — Explicit a-\i Scheme (Ref. [6]) Result 



53 



I.Or 


computed u 
scaled error 
exact solution 



-m.oio 


H 0.000 


H -0.010 


0.0 0.2 0.4 0.6 0.8 1.0 

Figure 1 8 — Dual-Mesh Explicit a-\i Scheme : 
Boundary Layer, Re = 100 


I.Or 

0.8 

0.6 

0.4 

0.2 


computed u 
scaled error 
exact solution 


i -i 1.5 


1.0 


4 

A 


* 

Ax 


A A A A AAAA A AA AA A -A - 




0 . 0 - 

0.0 0.2 0.4 0.6 0.8 1 .0 

Figure 19 — Explicit a-\i Scheme (Ref. [6]) : 
Boundary Layer, Re = 100 


0.5 


-\o.o 


■0.5 

x 


54 


Scaled Error 1 Scaled Error 



REPORT DOCUMENTATION PAGE 


Form Approved 
OMBNo. 0704-0188 


Public reporting burden tor this collection of information is estimated to average 1 hour per response. Including the time for reviewing instructions, searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 
collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson 
Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington. DC 20503. 


1. AGENCY USE ONLY (Leave btanfq 


4. TITLE AND SUBTITLE 


2. REPORT DATE 

November 1997 


3. REPORT TYPE AND DATES COVERED 

Technical Memorandum 


5. FUNDING NUMBERS 


The Implicit and Explicit a-/J. Schemes 


6. AUTHOR(S) 


WU-523-26-23-00 


Sin-Chung Chang and Ananda Himansu 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESSES) 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135-3191 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 


E-11005 


9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESSES) 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 


10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 


NASA TM— 97-206307 


11. SUPPLEMENTARY NOTES 


Sin-Chung Chang, NASA Lewis Research Center and Ananda Himansu, Cleveland Telecommunications Corporation, 
Solon, Ohio 44139. Responsible person, Sin-Chung Chang, organization code 5880, (216) 433-5874. 


12a. DISTRIBUTION/AVAILABILITY STATEMENT 


12b. DISTRIBUTION CODE 


Unclassified - Unlimited 
Subject Categories: 64 and 34 


Distribution: Nonstandard 


This publication is available from the NASA Center for Aerospace Information, (301) 621-0390. 


13. ABSTRACT (Maximum 200 wotxtm) 

Artificial numerical dissipation is an important issue in large Reynolds number computations. In such computations, the artificial dissipation 
inherent in traditional numerical schemes can overwhelm the physical dissipation and yield inaccurate results on meshes of practical size. In the 
present work, the space-time conservation element and solution element method is used to construct new and accurate numerical schemes such that 
artificial numerical dissipation will not overwhelm physical dissipation. Specifically, these schemes have the property that numerical dissipation 
vanishes when the physical viscosity goes to zero. These new schemes therefore accurately model the physical dissipation even when it is 
extremely small. The method of space-time conservation element and solution element, currently under development, is a nontraditional numerical 
method for solving conservation laws. The method is developed on the basis of local and global flux conservation in a space-time domain, in which 
space and time are treated in a unified manner. Explicit solvers for model and fluid dynamic conservation laws have previously been investigated. 
In this paper, we introduce a new concept in the design of implicit schemes, and use it to construct two highly accurate solvers for a 
convection-diffusion equation. The two schemes become identical in the pure convection case, and in the pure diffusion case. The implicit schemes 
are applicable over the whole Reynolds number range, from purely diffusive equations to purely inviscid (convective) equations. The stability and 
consistency of the schemes are analyzed, and some numerical results are presented. It is shown that, in the inviscid case, the new schemes become 
explicit and their amplification factors are identical to those of the Leapfrog scheme. On the other hand, in the pure diffusion case, their principal 
amplification factor becomes the amplification factor of the Crank-Nicolson scheme . We also construct an explicit solver with the treatment of 
diffusion being based on that in the implicit solvers. The explicit solver has only a CFL stability limitation on the Courant number, yet it retains the 
second-order spatial accuracy of the implicit schemes. 


14. SUBJECT TERMS 

Space-time; Conservation element; Solution element 1 16 . price code 


17. SECURITY CLASSIFICATION 
OF REPORT 

Unclassified 


18. SECURITY CLASSIFICATION 
OF THIS PAGE 

Unclassified 


19. SECURITY CLASSIFICATION 
OF ABSTRACT 

Unclassified 


20. LIMITATION OF ABSTRACT 


Standard Form 298 (Rev. 2-89) 

Prescribed by ANSI Std. Z39-18 
298-102 






















